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Abstract 

In  December  1996,  a  project  was  initiated  at  the  Institute  for  Systems  Research  (ISR),  under  an 
agreement  between  Northrop  Grumman  Electronic  Sensors  and  Systems  Division  (ESSD)  and  the  ISR, 
to  investigate  the  epitaxial  growth  of  silicon-germanium  (Si— Ge)  heterostructures  in  a  commercial  rapid 
thermal  chemical  vapor  deposition  (RTCVD)  reactor.  This  report  provides  a  detailed  account  of  the  ob¬ 
jectives  and  results  of  work  done  on  this  project  as  of  September  1997.  The  report  covers  two  main  topics 
-  modeling  and  model  reduction.  Physics-based  models  are  developed  for  thermal,  fluid,  and  chemical 
mechanisms  involved  in  epitaxial  growth.  Experimental  work  for  model  validation  and  determination  of 
growth  parameters  is  described.  Due  to  the  complexity  and  high  computational  demands  of  the  models, 
we  investigate  the  use  of  model  reduction  techniques  to  reduce  the  model  complexity,  leading  to  faster 
simulation  and  facilitating  the  use  of  standard  control  and  optimization  strategies.  Some  of  the  contents 
of  this  report  are  contained  in  [34] . 
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1  Introduction 


In  December  1996,  a  project  was  initiated  at  the  Institute  for  Systems  Research  (ISR),  under  an  agreement 
between  Northrop  Grumman  Electronic  Sensors  and  Systems  Division  (ESSD)  and  the  ISR,  to  investigate 
the  epitaxial  growth  of  silicon-germanium  (Si-Ge)  heterostructures  in  a  commercial  rapid  thermal  chemical 
vapor  deposition  (RTCVD)  reactor.  This  report  provides  a  detailed  account  of  the  objectives  and  results  of 
work  done  on  this  project  as  of  September  1997.  Some  of  the  contents  of  this  report  are  contained  in  [34]. 

1.1  Objectives 

Epitaxial  growth  of  Si-Ge  heterostructures  on  a  silicon  substrate  is  an  area  of  great  current  interest  [9,  15, 
17,  32].  This  is  mainly  due  to  the  superior  electrical  performance  and  manufacturing  economies  associated 
with  Si-Ge  devices  [9,  49]. 

One  such  epitaxial  growth  process  of  immediate  interest  to  Northrop  Grumman  ESSD  is  Si-Ge  co¬ 
deposition  on  a  patterned  silicon  substrate.  The  deposition  is  performed  using  a  commercial  RTCVD  reactor 
manufactured  by  Advanced  Semiconductor  Materials  (ASM),  Inc.  and  located  at  Northrop  Grumman  ESSD. 
The  patterned  substrate  in  this  case  is  a  silicon  wafer  upon  which  has  been  deposited  a  geometrically 
patterned  layer  of  silicon  oxide.  One  important  property  of  the  oxide  pattern  geometry  is  pattern  pitch , 
which  is  the  distance  between  peaks  in  the  oxide  layer  pattern.  Pattern  pitch  can  be  thought  of  as  a  measure 
of  how  densely  packed  the  pattern  geometry  is.  Its  importance  stems  from  the  depletion  phenomenon 
associated  with  tightly  packed  geometries  which  affects  the  uniformity  of  the  deposited  film. 

The  overall  objective  of  this  project  is  to  improve  the  manufacturing  effectiveness  of  the  epitaxial  growth 
process  and  the  RTCVD  reactor  by 

•  improved  understanding  of  the  processes  and  equipment  via  physical  and  mathematical  modeling,  and 

•  using  the  resulting  validated  models  for  optimization  of  process  conditions. 

More  specifically,  the  project  seeks  to  improve  product  quality  as  measured  by  deposition  uniformity. 
Currently,  the  process  engineer  operating  the  ASM  Epsilon-1  reactor  can  achieve  an  acceptable  level  of 
deposition  uniformity  (1%  variation)  on  bare  and  patterned  wafers,  but  the  uniformity  control  task  is  un¬ 
dertaken  on  mainly  a  trial-and-error  basis.  Furthermore,  the  limits  to  achieving  deposition  uniformity  in 
the  presence  of  microfeatures,  or  oxide  patterns,  on  the  wafer  surface  are  not  understood.  An  improved 
understanding  of  the  processes  and  equipment  will  provide  a  method  to  improve  or  optimize  process  settings 
and  conditions. 

Thus,  a  main  goal  is  to  develop  physics-based  models  for  thermal,  fluid,  and  chemical  mechanisms  in 
order  to  predict  deposition  characteristics  under  different  process  settings  and  conditions  such  as  pressure, 
flow  rate,  and  temperature.  The  modeling  effort  takes  into  account  lamp  characteristics,  reactor  geometry, 
chemical  reaction  kinetics,  and  optical  and  thermal  properties  of  the  materials  involved. 

Due  to  the  complexity  and  high  computational  demands  of  these  models,  we  are  investigating  the  use 
of  model  reduction  techniques  to  reduce  the  model  complexity,  leading  to  faster  simulation  and  facilitating 
the  use  of  standard  control  and  optimization  strategies.  An  important  objective  is  to  develop  new  model 
reduction  techniques  for  a  class  of  nonlinear  systems  that  are  exemplified  by  the  models  described  in  this 
report. 

1.2  Participants 

The  participants  in  this  project  are  as  follows:  the  ISR  team  consists  of  Andrew  Newman  and  Prof.  P.  S. 
Krishnaprasad,  and  the  Northrop  Grumman  team  consists  of  Sam  Ponczak,  Michael  O’Loughlin,  and  Paul 
Brabant.  The  ISR  participants  gratefully  acknowledge  the  assistance  of  the  Northrop  Grumman  participants 
and  the  use  of  their  equipment  for  experimental  purposes. 

1.3  Equipment 

The  equipment  used  for  thin  film  deposition  in  this  project  is  the  Epsilon-1  Reactor  System  manufactured 
by  ASM  Inc.,  Phoenix,  AZ.  One  of  these  reactor  systems  is  currently  used  as  a  production  tool  by  Northrop 
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Process  Chamber 


Figure  1:  ASM  Epsilon-1  RTCVD  reactor  wafer  handling  components.  Process  chamber  is  physical  area  of 
interest  for  deposition  process.  (Source:  ASM  Epsilon-1  Reactor  Manual.) 


Grumman  ESSD,  Linthicum,  MD.  The  equipment  and  process  models  described  in  this  report  are  specifically 
tailored  for  this  particular  reactor,  although  the  physical  principles  on  which  they  are  based  hold  for  similar 
reactors  and  mechanisms.  Also,  all  model  validation  and  other  experiments  performed  for  this  project  have 
taken  place  and  will  take  place  in  the  future  using  this  reactor. 

The  ASM  Epsilon-1  is  a  single-wafer,  lamp-heated,  RTCVD  reactor.  Figure  1  shows  a  schematic  of  the 
wafer  handling  components.  The  process  chamber  is  the  physical  area  of  interest  for  the  deposition  process. 
For  equipment  details  see  [5]. 

Figures  2  and  3  show  cross-sectional  views  of  the  reactor  section  and  the  lamp  assembly.  Process  gases 
flow  horizontally  through  the  process  chamber  from  inlet  slits  to  the  exhaust.  The  wafer  is  heated  by  an 
upper  and  lower  array  of  linear  tungsten-halogen  lamps  as  well  as  four  spot  lamps  directed  at  the  center 
of  the  wafer.  The  upper  lamp  array  illuminates  the  top  side  of  the  wafer  while  the  lower  array  and  the 
spot  lamps  illuminate  the  bottom  side  of  the  susceptor.  Four  thermocouples  at  the  center,  side,  front 
(upstream),  and  rear  (downstream)  of  the  susceptor  measure  temperature.  Thermocouple  readings  are  used 
in  proportional-integral-derivative  (PID)  control  loops  for  temperature  control.  The  reactor  can  operate  at 
both  atmospheric  and  reduced  pressure  (e.g.,  20  Torr).  The  graphite  susceptor  rotates,  typically  at  35  rpm. 
Radiation  from  the  wafer  edge  is  reduced  by  a  guard  ring.  The  chamber  walls  are  quartz  and  are  air  cooled. 

Northrop  Grumman  provided  two  tools  for  measuring  film  thickness  in  our  experiments,  nanospec  and 
ellipsometer.  The  Nanometrics  210  XP  Scanning  UV  Nanospec/DUV  Microspectrophotometer,  or  nanospec, 
is  an  instrument  for  measuring  the  thickness  of  optically  transparent  thin  (10  to  4000  nm)  films  on  silicon 
wafers.  The  speed  of  the  nanospec  makes  it  an  attractive  tool  for  taking  a  large  number  of  measurements. 
A  thickness  measurement  is  recorded  within  2  seconds  using  the  nanospec,  compared  with  approximately  15 
seconds  for  the  ellipsometer.  However,  the  nanospec  does  not  measure  poly  silicon  thickness  that  is  less  than 
100  Angstroms.  For  films  less  than  100  Anstroms,  we  used  the  ellipsometer. 

1.4  Project  Evolution 

RTCVD  of  Si-Ge  heterostructures  on  a  patterned  wafer  is  a  complicated  process  requiring  models  and 
analyses  that  focus  on  different  aspects  of  the  problem.  Therefore,  the  project  has  evolved  gradually  by  first 
considering  the  least  complex  aspects  and  tasks  and  then  building  upon  this  foundation  to  attack  the  more 
complicated  problems. 
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Upper  Susceptor  Exhaust 


Figure  2:  Cross-sectional  view  of  reactor  section  of  ASM  Epsilon-1  including  process  chamber.  Gases  flow 
from  inlet  injector  (left)  to  exhaust  (right).  (Source:  ASM  Epsilon-1  Reactor  Manual.) 


Figure  3:  Lamp  assembly  of  ASM  Epsilon-1.  Upper  array  illuminates  top  surface  of  wafer.  Bottom  array 
and  spot  lamps  illuminate  bottom  surface  of  susceptor.  (Source:  ASM  Epsilon-1  Reactor  Manual.) 
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Initial  efforts  in  this  project  focused  on  thermally  activated  epitaxial  growth  of  silicon  on  bare  silicon 
wafers.  Simulations  were  accomplished  by  hard  coding  stand-alone  programs  using  an  appropriate  program¬ 
ming  language  and  by  programming  commercially  available  general  purpose  and  specialized  computational 
fluid  dynamics  (CFD)  software  packages.  Both  of  these  software  implementation  approaches  numerically 
integrate  high-order  ordinary  differential  equation  (ODE)  approximations  of  partial  differential  equation 
(PDE)  models  for  the  relevant  balance  equations.  Models  that  incorporate  more  complicated  mechanisms 
associated  with  Si-Ge  co-deposition,  patterned  wafers,  and  deposition  in  the  mass  transport  controlled 
regime  are  under  investigation.  Different  software  implementations  of  models  are  being  integrated  to  work 
together  effectively.  Model  reduction  techniques  are  being  used  to  find  low-order  ODE  models  that  provide 
good  approximations  to  the  original  high-order  ODE  and  PDE  models. 

The  following  illustrates  a  summary  of  the  project  evolution. 

Si  Epitaxy  — >  Si-Ge  Epitaxy 

Bare  Wafers  — >  Patterned  Wafers 


Thermally  Activated 


Mass  Transport  Controlled 


Stand  Alone,  Hard  — »  Integrated  With  CFD 

Coded  Simulations 


PDE  and  High  — >  Low  Order  ODE  Models 

Order  ODE  Models 


1.5  Report  Outline 

This  report  is  organized  into  sections  as  follows.  Section  2  describes  in  detail  the  modeling  effort  including  the 
process  and  equipment  models  we  have  developed  and  experimental  work  for  validation  and  determination 
of  growth  parameters.  Section  3  discusses  our  work  in  model  reduction  as  it  is  applied  to  the  models  used  in 
this  project.  Section  4  outlines  future  directions  in  optimization  and  control,  and  presents  some  concluding 
remarks. 


2  Physics  Based  Models 

As  stated  in  Section  1,  one  goal  of  our  research  is  to  model  and  understand  the  phenomena  that  influence 
deposition  uniformity  in  the  presence  of  microfeatures  or  patterns  on  the  silicon  substrate.  The  first  task 
we  undertake  in  achieving  this  goal  is  the  development  of  physics-based  and  other  types  of  models  for 
processes  and  equipment.  These  models  will  provide  an  improved  understanding  of  the  process  and  equipment 
dynamics  and  a  foundation  for  subsequent  efforts  in  simulation,  optimization,  and  control. 

Development  of  sophisticated  models  will  progress  in  stages  of  increasing  generality,  fidelity,  and  complex¬ 
ity.  Ultimately,  the  models  must  retain  high  fidelity  for  a  wide  range  of  operating  conditions  for  temperature, 
pressure,  species  concentration,  and  flow  rate.  In  addition,  it  is  hoped  that  they  will  accurately  predict  de¬ 
position  characteristics  for  Si-Ge  thin  films  on  patterned  wafers.  Initially,  however,  we  make  simplifying 
assumptions,  restrict  models  to  only  the  most  important  processes,  and  consider  less  complicated  deposition 
tasks.  For  example,  we  begin  by  considering  deposition  of  thin  epitaxial  layers  of  silicon  on  a  bare  silicon 
substrate  under  the  assumption  that  kinetics  are  strongly  surface-reaction  controlled.  This  modeling  effort 
provides  a  foundation  for  a  series  of  minor  and  major  enhancements. 

2.1  Modeling  Overview 

The  general  modeling  scheme  is  illustrated  in  Figure  4  (see  [26]  for  a  similar  approach).  The  two  main 
components  of  the  modeling  effort  are  a  macroscopic  level  process  and  equipment  model,  and  a  microscopic 
level  feature  scale  model.  These  two  components  cooperate  to  achieve  the  overall  modeling  goals,  but 
separately  they  accomplish  qualitatively  different  tasks.  Therefore,  they  have  different  inputs,  outputs, 
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Figure  4:  Overview  of  modeling  effort  for  epitaxial  growth  in  ASM  Epsilon-1  reactor. 


implementation  strategies,  and  underlying  physics.  Cooperation  is  accomplished  by  using  certain  outputs  of 
the  macroscopic  level  model  as  inputs  to  the  microscopic  level  model. 

The  process  and  equipment  model  describes  phenomena  that  occur  at  a  macroscopic  level.  This  refers 
to  continuum  transport  of  energy,  mass,  and  momentum  in  the  solids  and  gases  that  comprise  the  reactor 
and  process  materials,  and  chemical  mechanisms  viewed  in  the  aggregate  sense.  The  motion  or  deposition 
of  individual  atoms  is  not  considered. 

The  feature  scale  model  describes  phenomena  that  occur  at  a  microscopic  level.  We  are  interested  in 
deposition  of  atoms  and  formation  of  alloys  on  a  patterned  substrate.  Therefore,  the  microscopic  model 
predicts  atom  motion  on  the  wafer  surface,  and  in  the  boundary  layer  just  above  the  surface,  under  given 
process  conditions,  as  provided  by  the  macroscopic  level  model. 

2.2  Modeling  Issues 

Epitaxial  growth  of  Si-Ge  thin  films  in  the  ASM  Epsilon-1  reactor  is  a  complicated  process  due  to  several 
factors.  The  deposition  kinetics  are  strongly  dependent  on  the  material  to  be  deposited,  source  gases, 
operating  conditions,  reactor  design,  and  wafer  properties.  Reactants  are  transported  to  the  wafer  via  gas 
flows  which  may  exhibit  swirling  or  turbulent  behavior  under  certain  conditions.  Heat  energy  is  exchanged 
among  lamps,  wafer,  chamber  walls,  and  flowing  gases  in  a  complex  manner.  Properties  of  the  materials 
involved  depend  on  process  variables  that  are  changing,  often  in  complicated  ways.  We  state  here  some 
complicating  issues  which  may  be  important  to  consider  for  development  of  high-fidelity  models,  and  some 
simplifying  assumptions  we  can  reasonably  adopt  for  dealing  with  them. 

2.2.1  Complications 

1.  The  wafer  rests  on  a  rotating  susceptor.  This  rotation  gives  rise  to  gas  flow  vorticity  and  dynamics  in 
the  axial  direction.  This  phenomenon  together  with  asymmetries  in  the  reactor  geometry  may  cause 
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azimuthal  symmetry  of  flow  over  the  wafer  surface  to  disappear.  Thus,  a  full  three-dimensional  model 
may  be  needed  to  predict  the  gas  flow  field. 

2.  The  ASM  Epsilon-1  reactor  can  operate  at  both  atmospheric  pressure  (AP)  and  reduced  pressure 
(RP).  Choice  of  pressure  mode  will  affect  deposition  kinetics.  Furthermore,  flow  can  be  turbulent  in 
the  AP  regime  [51],  although  it  has  been  observed  to  always  be  laminar  in  the  ASM  reactor  [29]. 

3.  Transport  of  reactant  species  to  the  wafer  surface  occurs  via  multi-component  diffusion  and  convection. 
Hence,  complicated  gas  flows  and  boundary  layer  effects  enter  into  the  deposition  kinetics  equations 
[51].  Depletion  effects  (reduction  in  reactant  species  concentration  at  the  downstream  end  of  the  flow) 
may  also  occur. 

4.  Deposition  rates  are  affected  by  germanium  content  in  the  reactant  gases.  This  phenomenon  is  the 
subject  of  much  recent  research  [22,  23,  25,  31],  which  provides  empirical  data  and  rules  of  thumb  but 
no  explicit  mathematical  models.  Effects  of  multiple  reactant  gases  complicate  the  deposition  kinetics 
models. 

5.  Reaction  kinetics  are  sensitive  to  patterns  (microfeatures)  on  the  wafer  surface.  For  example,  it  has 
been  shown  experimentally  that  thickness  dependence  on  an  oxide  pattern  is  a  strong  function  of 
pressure  [21]  and  deposition  rates  and  germanium  fraction  are  affected  by  pattern  pitch  and  closeness 
to  pattern  edges  [23]. 

6.  The  type  of  growth  (epitaxial,  poly  crystalline,  or  amorphous)  is  sensitive  to  the  operating  conditions 
[20,  30,  39,  44,  51],  most  importantly  temperature,  pressure,  flow  rate,  and  type  and  concentration 
of  reactant  species.  It  is  also  strongly  influenced  by  the  properties  of  the  substrate  surface  [29]. 
The  literature  provides  experimental  results  indicating  some  combinations  of  conditions  that  result  in 
epitaxial  growth  but  little  in  the  way  of  mathematical  models. 

7.  Deposition  can  be  either  mass  transport  or  surface  reaction  limited,  depending  on  wafer  surface  tem¬ 
perature.  The  transition  temperature  from  one  regime  to  the  next  is  a  function  of  silicon  precursor.  In 
the  ASM  reactor,  it  has  been  observed  that  for  silane,  mass  transport  limited  deposition  occurs  above 
850  C,  but  the  transition  temperature  is  higher  for  dichlorosilane,  tricholorsilane,  and  silicontetrachlo- 
ride  [29]. 

8.  Growth  rate  parameters  often  need  to  be  determined  experimentally.  Moreover,  it  is  likely  that  the 
measured  parameters  will  only  be  useful  for  deposition  of  one  particular  material  or  a  limited  range  of 
operating  conditions.  Ideal  Arrhenius  laws  must  be  modified  to  account  for  leveling  at  high  temperature 
and  for  pressure  dependence. 

9.  The  wafer  is  heated  with  upper  and  lower  lamp  arrays  which  are  divided  into  ten  lamp  groups  and  four 
lamp  heating  zones,  each  driven  with  separate  actuation.  Finding  the  relationship  between  a  particular 
lamp  group  or  zone  power  setting  and  radiant  heat  flux  intensity  on  the  wafer  surface  requires  either 
an  empirically  determined  heat  flux  intensity  profiles  (using  experimental  data  and  further  analysis), 
a  ray-trace  algorithm,  or  a  view  factor  model  (see,  e.g.,  [16,  24]).  Any  analytical  approach  needs  to  be 
verified  experimentally. 

10.  Heat  transfer  to  and  from  the  wafer  includes  conductive,  convective,  and  radiative  transport  mecha¬ 
nisms.  Heat  energy  is  exchanged  among  lamps,  wafer,  chamber  walls,  reflectors,  and  flowing  gas  in 
a  complex  manner.  High-fidelity  models  require  incorporation  of  chamber  wall  geometry  and  ma¬ 
terial  properties  for  computation  of  radiative  terms  due  to  reflections  and  nonuniformity  in  ambient 
temperature  [50]. 

11.  Material  properties  are  functions  of  the  process  variables,  whose  time  evolution  in  turn  depend  on  the 
material  properties.  For  example,  heat  transfer  in  the  wafer  depends  on  the  mass  density,  heat  capacity, 
thermal  conductivity,  and  emissivity  of  the  wafer.  Meanwhile,  mass  density  and  thermal  conductivity 
are  nonlinear  functions  of  the  wafer  temperature,  and  emissivity  is  a  nonlinear  function  of  deposition 
thickness.  Imperfections  in  the  wafer  may  cause  spatial  variation  in  properties.  Furthermore,  the  wafer 
itself  has  a  finite  thickness  so  that  there  may  exist  heat  transfer  in  the  axial  direction. 
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12.  Thermocouple  sensors  are  not  in  contact  with  the  wafer,  but  instead  are  located  at  four  locations  on 
the  susceptor  ring  and  susceptor  center.  If  temperature  measurements  are  used  in  a  model,  this  must 
be  taken  into  account. 

2.2.2  Simplifications 

Homoepitaxy  of  silicon  thin  films  on  a  silicon  substrate  is  a  less  complicated  and  better  understood  process 
than  heteroepitaxy  of  Si-Ge  thin  films.  Therefore,  the  initial  modeling  effort  focuses  on  homoepitaxy  of 
silicon  on  a  silicon  substrate  which  provides  a  foundation  for  a  progression  of  more  accurate  and  complex 
models.  Thus,  initially  we  eliminate  the  need  to  consider  effects  of  germanium  and  its  source  gases  on  the 
deposition  kinetics.  We  choose  silane  as  the  source  gas  for  homoepitaxy  of  silicon  since  this  provides  well 
known  and  simple  chemical  reactions,  with  reduction  occurring  at  lower  temperatures  than  for  the  other 
precursor  gases  such  as  silicon  tetrachloride  [51]. 

A  portion  of  the  research  described  in  this  report  focuses  on  thermally  activated  growth,  i.e.,  deposition 
for  which  reaction  kinetics  are  strongly  surface-reaction  controlled.  This  assumption  will  not  be  valid  when 
high  operating  temperatures  are  used.  However,  Northrop  Grumman  ESSD  performs  some  low  temperature 
epitaxial  growth  of  undoped  silicon  at  temperatures  between  600  C  and  800  C  [38].  In  that  portion  of 
our  work  where  this  assumption  is  invoked,  we  can  essentially  ignore  the  gas  flow  dynamics  which  control 
convective  transport  of  species  to  the  wafer  surface.  Instead,  we  can  assume  a  fixed  concentration  profile  of 
reactant  species  at  the  wafer  surface  as  a  parameter  in  the  kinetics  equations. 

Furthermore,  ignoring  gas  flow  dynamics  also  has  the  effect  that  wafer  rotation  is  no  longer  a  complicating 
factor.  In  fact,  wafer  rotation  helps  to  simplify  matters  by  eliminating  any  azimuthal  asymmetry  in  radiant 
heat  transfer  from  the  lamp  banks  to  the  wafer.  Thus,  the  wafer  thermal  dynamics  can  be  formulated  solely 
in  the  radial  direction.  Moreover,  since  a  susceptor  ring  is  used  we  can  realistically  ignore  wafer  edge  effects, 
i.e.,  lamp  heat  flux  incident  upon,  or  additional  radiative  losses  from,  the  side  of  the  wafer  edge. 

Initial  models  will  be  for  a  perfectly  flat,  perfectly  cylindrical,  homogeneous,  unpatterned,  silicon  wafer 
of  negligible  thickness.  Effects  of  microfeatures  will  be  considered  later.  Physical  constants  such  as  mass 
density,  heat  capacity,  thermal  conductivity,  and  emissivity  are  initially  assumed  constant,  i.e.,  no  variation 
with  temperature,  film  thickness,  position,  or  time. 

We  employ  a  view  factor  model  for  lamp  heating.  Radiative  heat  energy  exchange  within  the  chamber 
is  limited  to  only  that  between  the  wafer  surface  and  a  uniform  ambient,  in  this  case  the  chamber  walls. 
Reflections  off  chamber  walls  and  temperature  nonuniformity  throughout  the  chamber  and  walls  is  not 
included. 

The  physical  mechanisms  of  other  actuators  and  sensors  such  as  mass  flow  controllers  and  thermocouples 
are  not  modeled.  Instead  it  is  simply  assumed  that  they  function  as  intended  without  error  or  corruption. 

2.3  Process  and  Equipment  Models 

The  macroscopic  level  process  and  equipment  model  predicts  time  evolution  for 

•  heat  transfer  within  and  among  the  wafer,  chamber  walls,  and  process  gases  (temperature  fields  in 
solids  and  gases), 

•  momentum  transport  in  the  gas  phase  (gas  flow  velocity  vectors), 

•  mass  transport  in  the  gas  phase  (species  concentrations),  and 

•  chemical  reaction  kinetics  in  the  gas  phase  and  at  the  wafer  surface  (reaction  rates,  deposition  thick¬ 
ness). 

Thus,  it  consists  mainly  of  balance  equations  for  conservation  of  energy  (e.g.,  heat  equation),  and  momentum 
and  mass  (e.g.,  Navier-Stokes),  along  with  equations  that  describe  the  relevant  chemical  mechanisms  (e.g., 
Arrhenius  laws).  Boundary  conditions,  source  terms,  required  parameters,  and  any  other  details  of  the  model 
are  determined  using  available  data  regarding  the  reactor  geometry,  lamp  characteristics,  properties  of  the 
gases  and  solids  involved,  and  chemical  mechanisms.  In  the  course  of  this  project  some  of  this  data  has 
been  experimentally  measured  by  the  participants,  some  has  been  analytically  determined  via  mathematical 
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modeling,  and  some  has  been  gathered  from  appropriate  references  and  databases.  Each  particular  source 
of  data  will  be  indicated  in  this  report  in  the  appropriate  section. 

Inputs  to  the  model  consist  of  pressure,  inlet  flow  rates,  inlet  injector  slit  widths,  choice  of  gases,  and 
lamp  power  settings.  For  thermally  activated  growth  (surface  reaction  controlled),  the  heat  lamps  provide 
the  control  actuation  via  radiative  heat  transfer  to  the  wafer.  For  mass  transport  controlled  growth,  control 
actuation  is  achieved  by  adjusting  boundary  conditions  at  the  inlet. 

The  measured  outputs  of  the  model  include  aggregate  growth  rate,  spatial  uniformity,  and  film  composi¬ 
tion.  These  are  the  variables  which  we  ultimately  wish  to  control.  We  also  assume,  when  convenient,  access 
to  other  variables,  such  as  the  temperature  at  given  points  on  the  wafer  surface,  which  is  used  an  an  input 
to  the  feature  scale  model.  These  variables  may  or  may  not  be  available  as  measured  outputs.  If  in-situ 
measurement  is  not  possible  then  for  real-time  control  a  method  of  estimating  them  would  be  required. 

2.3.1  Implementation 

As  stated  in  the  previous  section,  the  macroscopic  level  process  and  equipment  model  consists  of  evolution 
equations  and  associated  boundary  conditions  for  transport  of  heat,  mass,  and  momentum  along  with  gas 
phase  and  surface  chemical  reaction  kinetics.  We  take  two  approaches  to  implementation.  One  approach  is  to 
hard-code  individual  stand-alone  software  modules  in  a  suitable  programming  language  such  as  MATLAB1 
to  simulate  individual  components  of  the  model.  These  individual  components  can  then  be  interfaced  to 
simulate  larger  portions  of  the  model  or  the  overall  model.  The  other  approach  is  to  use  commerically  avail¬ 
able  general  purpose  computational  fluid  dynamics  (CFD)  software  packages  with  capability  for  simulation 
of  complex  fluid  flow,  heat  transfer,  and  chemical  reaction  systems.  In  addition,  supplemental  code  can  be 
developed  for  integration  with  the  CFD  software  for  modeling  phenomena  specific  to  the  equipment  and 
processes  of  interest.  We  shall  see  that  both  (hard-coding  and  CFD)  approaches  are  useful  and  necessary, 
and  that  the  overall  implementation  requires  some  degree  of  redundancy  between  the  two. 

There  currently  exist  commercially  available  general  purpose  CFD  software  packages  that  provide  the 
required  capabilities  for  implementing  the  model  including 

•  body  fitted  grids  for  discretization  of  complicated  physical  domains  such  as  those  found  in  a  CVD 
reactor, 

•  efficient  algorithms  for  numerically  integrating  the  fully  coupled  nonlinear  versions  of  the  transport 
equations,  and 

•  Monte-Carlo  ray-tracing  and  view  factor  algorithms  for  modeling  radiative  heat  transfer  among  sur¬ 
faces. 

We  are  using  two  of  these  general  purpose  CFD  packages,  Fluent2  and  PHOENICS3,  to  provide  high-fidelity 
simulations  of  the  fluid  flow,  heat  transfer,  species  transport,  and  chemical  mechanisms  in  the  ASM  Epsilon-1 
reactor.  One  attractive  feature  of  PHOENICS  is  a  specialized  add-on  program  called  PHOENICS-CVD  for 
modeling  CVD  processes  which  includes  specialized  databases  for  optical  and  thermal  properties  of  materials 
and  gas  phase  and  surface  chemistry  of  many  commonly  used  reactions.  It  is  not  our  purpose  here  to  provide 
a  complete  review  of  the  capabilities  of  these  software  packages.  We  will  indicate  when  and  how  they  are 
used  in  this  modeling  effort  as  the  need  arises  in  this  report. 

The  CFD  approach  does  not  do  everything  we  need.  Not  only  are  general  purpose  CFD  programs 
restricted  to  modeling  at  the  macroscopic  level,  but  the  output  data  is  not  particularly  suitable  for  integration 
into  a  model  at  the  microscopic  level.  For  patterned  wafers,  we  are  studying  mechanisms  at  the  feature  scale. 
We  need  models  which  focus  attention  on  what  is  happening  at  the  wafer  surface  and  which  provide  very 
high  spatial  resolution  there.  General  purpose  CFD  simulations  with  a  spatial  resolution  (discretization) 
suitably  fine  for  feature  scale  modeling  would  be  computationally  prohibitive. 

Furthermore,  implementation  of  feedback  control  and  optimization  techniques  with  the  reactor  model  as 
the  plant  would  be  cumbersome  and  impractical  using  the  general  purpose  CFD  model.  Model  reduction 

lrThe  Mathworks,  Inc.,  Natick,  MA 
2 Fluent,  Inc.,  Lebanon,  NH 
3 CHAM  Ltd.,  London,  UK 
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Figure  5:  Hard-coded  individual  MATLAB  programs  are  integrated  with  general  purpose  and  specialized 
CFD  models  to  produce  the  desired  overall  input-output  model  of  epitaxial  growth  the  ASM  Epsilon-1 
reactor. 

strategies  for  development  of  low-order  models  can  be  applied  to  the  hard-coded  stand-alone  programs  to 
provide  useful  approximations  of  the  original  models  for  use  in  optimization  and  control. 

There  is  a  necessary  degree  of  redundancy  between  the  two  approaches.  For  example,  heat  transfer 
in  the  solid  wafer  is  implemented  both  via  general  purpose  CFD  software  and  via  stand-alone  MATLAB 
programs.  In  the  CFD  approach  the  spatial  discretization  of  the  solid  wafer  is  relatively  coarse  but  suitable  for 
simulating  the  overall  heat  transfer  in  the  reactor.  In  the  MATLAB  implementation  the  spatial  discretization 
can  be  set  arbitrarily  fine  without  seriously  degrading  computational  performance,  and  the  resulting  code  is 
immediately  available  for  use  in  model  reduction,  optimization,  and  control  schemes. 

In  summary,  the  CFD  approach  yields  more  realistic  solutions,  since  the  sophisticated  software  packages 
solve  the  fully  coupled  nonlinear  versions  of  the  transport  equations.  It  would  be  impractical  for  the  authors 
to  attempt  to  replicate  all  of  these  advanced  features  into  our  stand-alone  hard-coded  programs.  Instead,  the 
hard-coded  programs  invoke  many  simplifying  assumptions  to  make  their  implementation  practical.  Both 
approaches  are  integrated  to  accomplish  the  project  objectives.  Figure  5  illustrates  a  strategy  for  integrating 
the  two  approaches  to  build  a  complete  input-output  model  for  describing  morphology  evolution  for  epitaxial 
growth  in  the  ASM  Epsilon-1  reactor. 

2.3.2  General  Purpose  and  Specialized  CFD 

As  stated  earlier,  we  are  using  two  general  purpose  CFD  software  packages,  Fluent  and  PHOENICS,  to 
model  fluid  flow,  heat  transfer,  and  chemical  reactions  in  the  ASM  Epsilon-1  reactor.  These  packages 
incorporate  up-to-date  modeling  techniques  and  a  wide  range  of  physical  models.  In  addition,  we  are  using 
a  specialized  add-on  program,  PHOENICS-CVD,  which  has  features  and  databases  of  special  use  for  CVD 
applications.  We  cannot  provide  a  full  discussion  of  the  equations,  algorithms,  discretization  techniques,  and 
other  characteristics  of  these  software  packages.  See  [11]  for  a  complete  description  of  Fluent  and  various 
articles  in  [48]  for  a  detailed  description  of  PHOENICS-CVD  and  several  example  applications. 

Here  we  demonstrate  the  CFD  modeling  approach  as  applied  to  the  ASM  Epsilon-1  reactor  by  presenting 
some  results  from  steady-state  simulations  of  certain  phenomena  of  interest.  These  simulations  assumed  the 
following  constant  operating  conditions:  wafer  temperature  725  C,  chamber  wall  temperature  425  C,  inlet 
gas  temperature  25  C,  and  inlet  flow  rate  1.5  slm  of  2%  silane  in  hydrogen.  Effects  due  to  surface  reactions 
are  included. 

A  2-dimensional  model  of  the  ASM  Epsilon-1  reactor  was  developed  using  a  Cartesian  geometry  and 
nonuniform  grid.  This  model  represents  a  2-dimensional  slice  of  the  process  chamber  parallel  to  the  direction 
of  flow,  i.e.,  from  inlet  to  outlet.  Symmetry  is  assumed  in  the  direction  perpendicular  to  the  plane  of  the 
slice. 

Figure  6  shows  the  steady-state  condition  for  gas  flow  velocity.  The  gas  velocity  magnitude  is  large  near 
the  inlet,  with  some  recirculation  vortices.  The  velocity  magnitude  is  small  as  it  passes  the  wafer,  with  the 
expected  parabolic  flow  profile.  Figure  7  shows  the  steady-state  condition  for  gas  temperature.  The  gas  is 
cold  as  it  leaves  the  inlet  and  is  heated  as  it  passes  the  wafer.  Figure  8  shows  the  steady-state  condition  for 
silane  mole  fraction.  Surface  reactions  cause  a  depletion  effect  as  the  concentration  of  silane  decreases  from 
upstream  to  downstream  side  of  the  wafer.  Figure  9  shows  the  mass  flux  of  deposited  silicon.  Deposition  is 
a  maximum  at  the  leading  edge  of  the  wafer  due  to  the  depletion  effect. 
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Figure  6:  Steady-state  gas  velocity  vectors  (fluid  flow)  in  2-dimensional  slice  of  ASM  Epsilon-1  process 
chamber. 
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Figure  7:  Steady-state  gas  temperature  in  2-dimensional  slice  of  ASM  Epsilon-1  process  chamber. 
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Figure  8:  Steady-state  silane  mole  fraction  in  2-dimensional  slice  of  ASM  Epsilon-1  process  chamber. 
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Figure  9:  Plot  of  deposited  silicon  mass  flux  versus  wafer  radial  position  for  ASM  Epsilon-1. 
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Figure  10:  Wafer  with  coordinates  and  measurements  used  in  the  heat  transfer  model. 


2.3.3  Wafer  Heat  Transfer 

The  model  for  wafer  heat  transfer  is  a  modified  version  of  models  presented  in  [2,  8,  28,  41,  42].  It  is  based 
on  an  energy  balance  for  a  heat  conducting  solid  which  emits  and  absorbs  heat  radiation  at  its  boundary 
surfaces.  The  model  takes  into  account  simplified  effects  of  conductive,  radiative  (including  lamp  heating), 
and  convective  heat  transfer.  Both  a  continuum  model  and  a  discretized  version  are  presented  here,  based 
on  identical  principles  of  energy  balance.  Although  the  wafer  is  a  continuous  solid  body,  a  discretized  model 
is  required  for  purposes  of  numerical  solution. 

Since  the  wafer  shape  is  assumed  to  be  a  perfect  cylinder,  the  model  is  formulated  in  cylindrical  coor¬ 
dinates  with  radial  variable  r,  azimuthal  variable  #,  and  axial  variable  z.  The  wafer  radius  is  denoted  Rw 
and  the  wafer  thickness  is  denoted  A2,  so  that  the  top  surface  of  the  wafer  has  ^-coordinate  Az  and  the 
bottom  surface  has  ^-coordinate  0.  Figure  10  shows  a  diagram  of  the  wafer  with  coordinates.  Note  that  for 
purposes  of  the  initial  analysis,  the  wafer  and  susceptor  have  been  combined  into  a  single  homogeneous  solid 
body. 

Continuum  Model 


The  temperature  field  in  the  solid  wafer  is  denoted  Tw  =  Tw(t,  r,  #,  z)  where  t  is  the  time  variable.  Time 
evolution  of  Tw  is  described  by  a  partial  differential  equation  (PDE)  (usually  referred  to  as  the  heat  equation) 
which  models  heat  conduction  within  the  wafer,  together  with  boundary  conditions  (BCs)  which  model  net 
heat  flow  to  and  from  the  wafer  boundary  surfaces  (top,  bottom,  and  edge).  The  PDE  is  given  in  cylindrical 
coordinates  by 


P  w  Cp 


dTw 

dt 


l±(k  rWA  +  ±±(k  \  +  JL(k 

r  dr  \  w  dr  )  r2  d9  \  w  dO  )  dz  \  w  dz  ) 


(1) 


for  t  >  0,  0  <  r  <  Rw ,  0  <  0  <  2w,  and  0  <  2  <  A,  where  pw  is  the  mass  density  of  the  wafer,  CPw  is  the 
heat  capacity  of  the  wafer  (the  product  Mw  =  pw  CPw  is  often  referred  to  as  the  wafer  thermal  mass) ,  and 
kw  is  the  thermal  conductivity  of  the  wafer.  The  associated  BCs  are  given  by 


dTw 


kyj 

kyj 

kyj 


dr 

dTw 

dr 

dTw 

dz 

dTw 

dz 


0 ;  r  =  0 

qedge(0,z);  T  =  Rw 
-qbottom(r,6);  2  =  0 

qtop(r,  0);  z  =  Az 


(2) 

(3) 

(4) 

(5) 


where  the  first  BC  results  from  symmetry  about  the  wafer  center,  and  qedge ,  Qbottom ,  and  qtop  represent 
the  net  heat  flow  per  unit  surface  area  to  and  from  the  wafer  edge,  bottom,  and  top  boundary  surfaces, 
respectively,  and  will  be  described  in  more  detail  later. 

For  purposes  of  modeling  film  growth,  we  focus  our  attention  on  the  top  surface  of  the  wafer  where 
reactions  take  place.  Invoking  the  assumption  of  azimuthal  symmetry,  so  that  no  temperature  gradients 
exist  in  the  azimuthal  direction  (i.e.,  8Tw/d0  =  0),  time  evolution  of  Tw  at  the  wafer  top  surface  is  described 
by 


Pw  Cj 


Pw 


dU, 

dt 


I A  (k 

r  dr  \  w  dr  J 


dTw\ 
dz  J 


(6) 
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for  t  >  0,  0  <  r  <  Rw ,  and  z  =  Az  with  BCs  remaining  the  same  as  above  except  qedge  =  Qedge(z )> 
q_bottom  =  Qbottomir ),  and  qtop  =  qtop(r).  We  also  assume  that  the  wafer  thickness  is  sufficiently  small  so  that 
no  thermal  gradients  exist  in  the  axial  direction  within  the  wafer  interior.  Therefore,  we  approximate  the 
axial  gradient  term  at  the  top  surface  by 


9_ 

dz 


1 

a! 


9TW  ,  7  dTw  , 

’  dz  U=A‘  w  dz  U=0 


( Qtop  Q bottom ) 


where  we  have  made  substitutions  using  the  appropriate  BCs.  The  resulting  PDE  describes  the  evolution  of 
the  wafer  top  surface  temperature  field  as  a  function  of  time  and  radial  position, 


with  BCs 


Pw  Cj 


Pw 


dTw 
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I A  (k  rdJA 

r  dr  \  w  dr  ) 


T  ^  ( qtop  qbottom ) 


9TW 
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dr 

dTw 

dr 


0  ;  r  =  0 

qedge(Az )  ;  r  =  Rw 


(7) 

(8) 

(9) 

(10) 


Now,  we  must  find  expressions  for  qtop  and  qbottom ,  the  net  heat  flow  into  the  top  and  bottom  surfaces  of 
the  wafer.  For  the  initial  analysis,  we  assume  that  the  top  and  bottom  surfaces  are  subject  to  identical  heat 
transfer  mechanisms,  and  let 

qtop  +  qbottom  —  q  +q  +q  +q  (iij 

where  the  terms  on  the  right  hand  side  represent  the  flow  of  thermal  energy  to  and  from  the  wafer  and  are 
dependent  on  time,  position,  and  wafer  temperature.  In  particular,  qem  is  radiative  energy  emitted,  qab  is 
radiative  energy  absorbed,  qconv  denotes  energy  losses  due  to  convective  heat  transfer,  and  qdlst  is  energy 
transfer  due  to  unmodeled  effects  such  as  heat  generated  by  chemical  reactions.  In  what  follows  we  ignore 

qdist 

The  term  qem  represents  radiative  losses  from  the  wafer.  We  assume  qab  depends  on  radiant  heat  flux 
from  a  uniform  ambient,  in  this  case  the  chamber  walls,  and  radiant  heat  flux  from  the  lamps,  but  without 
reflections  or  other  effects.  The  individual  terms  are  given  by 

„em  _  c\  rpA 

q  —  Z€w&b  -L  w 

10 

qab  =  2  awab  Tc4  +  aw  ^  QiUi 

where  denotes  the  Boltzmann  constant,  ew  denotes  the  wafer  emissivity,  aw  denotes  the  wafer  absorptivity, 
Tc  denotes  the  uniform  ambient  temperature  of  the  chamber  walls,  Qi  =  Qi(r)  is  a  function  of  position 
describing  the  heat  flux  intensity  incident  on  the  wafer  due  to  the  i-th  lamp  group,  and  Ui  =  Ui(t)  is  the 
time- varying  actuated  power  level  of  the  i-th  lamp  group. 

The  convective  term  is  given  by 

qconv  =  -hv(Tw-Tg)  (14) 

where  hv  denotes  the  convective  heat  transfer  coefficient  and  Tg  denotes  the  temperature  of  the  gas  flowing 
past  the  wafer.  Note  that  we  have  assumed  a  constant  uniform  gas  temperature.  In  order  to  estimate  hVl  we 
assume  flow  in  the  process  chamber  is  a  laminar  flow  along  a  flat  plate.  The  mean  heat  transfer  coefficient 
is  given  in  [37]  pp.  233-235  as 

hv  =  2  [0.332  kg  Fr1/3  (. Re 1/2 /L)}  (15) 

where  kg  denotes  the  gas  thermal  conductivity,  Pr  denotes  the  gas  Prandtl  number,  Re  denotes  the  gas 
Reynolds  number,  and  L  denotes  the  length  of  the  chamber.  We  have  computed  the  Reynolds  number  Re  to 


(12) 

(13) 
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be  approximately  27  for  the  flow  in  the  ASM  Epsilon-1  during  a  typical  deposition  run,  thus  confirming  the 
laminar  assumption.  The  calculated  value  of  hv  was  then  validated  using  flow  and  temperature  data  from  a 
corresponding  CFD  simulation.  See  Appendix  A  for  values  of  all  of  these  physical  constants. 

It  is  sometimes  convenient  to  assume  that  the  wafer  is  a  graybody,  so  that  ew  =  aw  for  all  relevant 
wavelengths  of  radiation  and  wafer  temperatures.  However,  we  do  not  make  this  assumption  here,  and  use 
different  values  for  emissivity  and  absorptivity.  We  also  note  that  the  parameters  pw ,  CPw ,  and  kw  can  be 
modeled  as  nonlinear  functions  of  TWl  and  the  parameters  ew  and  aw  can  be  modeled  as  nonlinear  functions 
of  Tw  and  deposition  thickness.  However,  we  invoke  the  assumption  that  mass  density,  heat  capacity,  thermal 
conductivity,  emissivity,  and  absorptivity  are  constant,  i.e.,  no  variation  with  temperature,  film  thickness, 
position,  or  time.  The  PDE  model  specializes  to 


dTw  =  kw  Id  f  dTw\ 
dt  pwCPw  r  dr  V  dr  ) 


PwCpw  A; 


<Tg 


Tw)  + 


2  (Jf)  €yj  /^4  rri4: 

PwCPw  A?  c  "  - 


OLw 


Pw  CPw  A 


10 

i=  1 


(16) 


where  we  recall  that  Tw  =  Tw(t,  r),  Qi  =  Qi(r),  and  Ui  =  Ui(t). 

Since  the  guard  ring  insulates  the  wafer  from  radiation  directed  at  its  edge  boundary  surface,  we  assume 
zero  heat  transfer  at  the  wafer  edge  so  that 

Qedge  0 

giving  the  boundary  conditions  (BCs) 


dTw 

dr 

9TW 

dr 


o 

II 
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Figure  11:  Wafer  discretized  into  annular  regions.  Heat  transfer  in  wafer  described  by  energy  balance  for 
each  discrete  element. 

Discretized  Model 

The  continuum  model  as  given  by  the  above  PDE  and  BCs  can  be  discretized  using  a  suitable  scheme,  e.g., 
finite  differences  or  finite  elements.  However,  for  our  simplified  model  it  is  easier  to  formulate  a  discretization 
by  applying  the  energy  balance  principles  directly  to  individual  annular  elements  of  the  wafer.  The  general 
idea,  which  divides  the  wafer  into  annular  regions,  is  illustrated  in  Figure  11.  Annular  regions  will  be 
numbered  from  1  to  n  with  element  1  being  the  innermost  disk  and  element  n  being  the  outermost  annular 
region.  Element  i  has  mean  radius  r(i),  and  is  bounded  by  an  outer  cylinder  of  radius  rout,  inner  cylinder 
of  radius  r^n,  top  surface  at  z  =  A2,  and  bottom  surface  at  z  =  0.  The  discretization  is  uniform  so  that 

Ar  =r(i)  —  r(j) 


is  constant  for  all  i,  j. 

The  usual  symmetry  assumptions  are  invoked  so  that  temperature  is  dependent  upon  radial  position  and 
time  only.  The  discretized  wafer  temperature  field  is  given  by  the  n-vector  Tw(t),  where  the  i-th  entry  of 
Tw(t)  represents  the  temperature  at  radial  position  r(i)  and  time  t. 

The  wafer  heat  transfer  model  is  then  given  by  the  ODE 

Tw  =  ACTW  +  ArT^  +  AVTW  +  T  +  BP  (19) 

where  Ac,  Ar,  and  Av  are  n  x  n  matrices  representing  the  effects  of  conductive,  radiative,  and  convective 
heat  transfer  mechanisms,  respectively,  T  is  a  constant  n-vector  which  accounts  for  the  gas  and  chamber 
wall  temperature,  B  is  a  n  x  m  matrix  of  discretized  lamp  zone  radiant  intensity  profiles,  and  P  =  P(t)  is 
a  m- vector  of  control  inputs  corresponding  to  lamp  zone  power  levels.  For  the  Epsilon-1  reactor,  there  are 
m  =  4  independently  actuated  lamp  zone  control  inputs.  We  present  the  details  of  the  ODE  model  below. 

The  top  surface  area,  volume,  and  mass  of  annular  region  i  are  given,  respectively,  by 


S(i) 

=  n  (rout(i)2  -  rin(i)2) 

V(i) 

=  S(i)  Az 

(20) 

m(i) 

=  pwV{i) 

The  matrix  representing  conductive  heat  transfer  is  then  represented  by  the  tridiagonal  matrix  given  by 
the  entries 


Ac(i,  i) 

Ac(i ,  i  +  1) 


Ac(i,  i  —  1) 


2  kw  ^out(^)  T  Tin(i) 

Pw  Cpw  Ar  r out  (i)2 

2  kw  T  out  (0 

Pw  Cpw  Ar  rout  (i)2 
2  kw  Tin  (i) 

Pw  Cpw  Ar  rout(i 


(21) 
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for  i  =  2, ...  n  —  1  and 


4  (1  ]\  —  2kw  1 

cl  ’  '  Pw  CPw  Ar  rout(l) 

4  (l  o')  _  2kw  1 

cl  ’  ^  Pw  CPw  Ar  rout(l) 

a  /  \  2  kw  Tin  (^) 

c[n,n)  ~  pw  CPw  Ar  rout(n)2  —  rin(n)2 

a  /  \  2  kw  Vin  (jl) 

Ac{n,n  1)  = 

pw  CPw  Ar  rout(n)2  -  rin(ny 

(22) 

where  we  note  that  zero  heat  flux  BCs  have  been  incorporated  into  the  model  via  boundary  elements  of 
matrix  Ac. 

The  matrices  representing  radiative  transfer  from  wafer  surface  to  chamber  walls  and  convective  heat 
transfer  from  the  process  gases  to  wafer  are  given,  respectively,  by 

Ar  =  diag[  ^  ] 

Pw  ^Pw  Z 

(23) 

and 

Av  =  diag[  ] 

Pw  ^pw  ^ Z 

(24) 

where  we  note  that  Ar  and  Av  take  diagonal  form  as  a  result  of  the  simplifications  made  in  our  model. 

The  effect  of  radiation  from  chamber  wall  to  wafer  and  convective  transfer  from  gas  to  wafer  are  incor¬ 
porated  into  the  constant  vectors  Tr  and  Tv  whose  entries  for  i  =  1 . . .  n  are  given  by 


r  r(i) 

r  „(*) 


6C  QLW  ^4 
pw  CPw  Az 


pw  CPw  A  2 


(25) 

(26) 


where  ec  is  the  emissivity  of  the  quartz  chamber  walls.  These  effects  are  combined  by  summing  into  one 
constant  vector 

r  =  rr  +  rv.  (27) 

Discretized  lamp  heat  flux  intensity  profiles  as  given  in  Section  2.3.4  are  arranged  in  a  matrix  Q  and 
incorporated  into  the  influence  matrix  B  given  by 


B  = 


Pw  Cpw  A; 


Q- 


(28) 


To  avoid  problems  of  scaling  in  computational  work,  we  normalize  variables  and  parameters  so  that 
all  units  cancel,  i.e.  write  the  model  in  dimensionless  form.  It  is  customary  to  adopt  a  notation  for  the 
dimensionless  variables,  e.g.  Tw  becomes  Tw.  Instead,  we  denote  the  dimensionless  variables  by  the  same 
symbol  as  their  dimensional  counterparts  and  caution  the  reader  to  keep  this  in  mind.  The  conversions  are 


Qi  ~ ^ 


Qi 

Qref 


r 

r  ~ *  ~R~' 

-Ll'W 


and  are  also  given  in  Appendices  B  and  C. 

It  has  been  observed  in  the  ASM  reactor  that  Tc,  the  chamber  wall  temperature,  is  approximately  300  K 
less  than  wafer  temperature  [29]  during  a  typical  processing  run.  As  reference  values  we  select  a  wafer 
temperature  of  1000  K  and  chamber  wall  temperature  of  700  K.  The  reference  thickness  href  of  1.0  micron 
was  selected  because  it  is  on  the  order  of  the  thickness  of  films  we  are  interested  in  growing.  The  reference 
heat  flux  Qref  of  29.24  W/cm2  was  computed  using  the  lamp  power  specification  of  6  kW  radiating  over 
one-half  of  a  spherical  surface  area  of  radius  2.25  inches.  All  values  are  given  in  Appendix  A. 
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Figure  12:  Organization  (top  view)  of  upper  and  lower  lamp  arrays  and  spot  lamps:  individual  lamps  are 
assigned  to  lamp  groups  and  heat  zones  as  shown. 


Using  the  dimensionless  variables  as  given  above,  the  parameters  (matrices  and  vectors)  in  equation  (19) 
become 


7fc  ~d2  ^  Zlr  -A  r  Tj?  Ar  ,  Av  -A  r  Av  ,  Tr  -A  —  Tr  ,  Tv  — >  — -  ,  B  -A  —  B 

Bw  ±c  J-c  2  c 

to  yield  an  ODE  model  equivalent  to  (19). 

The  ODE  (19)  can  be  numerically  integrated  to  approximate  the  PDE  (16),  given  appropriate  initial 
conditions,  TWo  =  Tw( 0),  to  determine  the  temperature  field  on  the  wafer  surface  as  a  function  of  time 
and  radial  position.  Typically,  the  initial  condition  is  a  constant  temperature  field  set  to  ambient,  i.e., 
Tw o  =  [700  ...  700] T.  For  the  work  described  in  this  report  we  used  a  fourth  and  fifth  order  Runga-Kutta 
integration  scheme  to  perform  the  numerical  integrations.  The  discretization  resolution  was  typically  set  at 
n  =  101. 

2.3.4  Lamp  Heating 

Heat  transfer  in  the  RTCVD  reactor  is  dominated  by  radiative  effects  such  as  radiative  transfer  between 
wafer  and  chamber  and  radiative  transfer  from  lamps  to  wafer.  The  importance  of  a  lamp  heating  model,  as 
one  component  of  the  overall  heat  transfer  model,  is  apparent,  especially  for  describing  a  thermally  activated 
deposition  process. 

Preliminaries 

Figure  12  illustrates  the  layout  and  organization  of  the  heat  lamps  in  the  ASM  Epsilon-1  reactor.  Upper 
and  lower  arrays  are  perpendicular  to  each  other,  with  spot  lamps  in  the  center.  Individual  lamps  are 
combined  into  groups,  which  are  further  combined  into  zones.  It  is  the  heat  zones  that  are  controlled 
independently  in  the  ASM  Epsilon-1.  Thus,  the  Epsilon-1  has  four  control  inputs  for  heating.  The  zone 
name  roughly  corresponds  to  the  area  of  the  wafer  that  receives  the  most  intense  illumination  from  the 
particular  zone,  e.g.  center,  front  (upstream),  rear  (downstream),  and  side.  The  number  of  lamps  in  each 
zone  varies  as  indicated. 

Lamp  heating  appears  in  the  dynamic  model  for  wafer  heat  transfer  as  the  control  input  term,  and  is 
described  by  a  set  of  influence  functions  which  represent  the  heating  effect,  or  incident  heat  flux  intensity, 
of  each  lamp  heat  zone  on  the  wafer  surface.  Each  influence  function  is  multiplied  by  a  corresponding  di¬ 
mensionless  scalar  control  input  which  represents  the  proportion  of  full  power  applied  to  the  particular  lamp 
zone  actuator.  Each  control  input  is  typically  a  function  of  time  and/or  sensor  measurements,  e.g.  thermo¬ 
couple  temperature  readings.  In  the  heat  transfer  model  (19)  the  lamp  influence  functions  are  appropriately 
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discretized  and  arranged  into  a  matrix,  denoted  B.  Similarly,  the  control  inputs  are  arranged  into  a  vector 
of  time-varying  functions,  denoted  P. 

The  influence  functions  physically  correspond  to  lamp  heat  flux  intensity  spatial  profiles.  A  heat  flux 
intensity  spatial  profile  is  a  function  of  position  on  the  wafer  surface  whose  value  at  a  spatial  point  is  the  heat 
flux  intensity  (measured  in  Watts  per  unit  area)  irradiated  on  the  given  point.  In  the  case  where  we  have 
azimuthal  symmetry,  such  as  when  the  wafer  rotates  at  a  uniform  rate,  azimuthal  variations  are  averaged 
over  360  degrees,  so  that  each  profile  is  a  function  of  radial  position  only.  Such  a  profile  is  determined  for  each 
individual  lamp  in  the  ASM  Epsilon-1  reactor.  Then,  individual  profiles  are  combined  using  superposition 
into  profiles  for  the  ten  lamp  groups  and  four  lamp  zones  of  the  Epsilon-1. 

The  flux  intensity  profiles  are  determined  analytically  and  verified  experimentally.  One  possible  experi¬ 
mental  approach  would  be  to  infer  flux  intensity  profiles  from  temperature  measurement  data.  For  example, 
an  instrumented  wafer  with  nine  attached  thermocouples  is  available  at  Northrop  Grumman  to  provide  tem¬ 
perature  data  at  nine  points  on  the  wafer  surface.  However,  based  on  past  experience,  there  is  some  doubt 
that  the  instrumented  wafer  can  provide  data  of  sufficient  accuracy  for  this  purpose,  especially  when  used 
under  typical  flow  conditions  in  the  reactor.  An  alternative  method  for  temperature  measurement  is  inferring 
temperature  from  deposition  thickness  data.  We  used  this  method  for  validation  purposes  and  it  is  described 
later.  The  main  problem  with  using  experimental  temperature  data  to  compute  lamp  heat  flux  intensity 
profiles  is  achieving  the  desired  spatial  resolution.  With  the  instrumented  wafer,  we  get  readings  at  only 
nine  spatial  points.  Using  the  deposition  thickness  approach,  a  very  large  number  of  thickness  measurements 
must  be  taken,  and  it  is  difficult  or  impossible  to  determine  the  exact  spatial  location  of  the  measurement  on 
the  wafer.  Therefore,  the  analytically  determined  profiles  are  used  in  the  model,  since  they  can  be  computed 
using  an  arbitrarily  fine  discretization  at  exactly  those  points  we  choose.  Experimental  results  are  used  only 
for  comparative  validation. 

In  what  follows  we  describe  our  methodology  for  determining  the  profiles,  present  some  preliminary 
results,  and  discuss  improvements  which  are  necessary  for  the  lamp  heating  models  to  accurately  predict  the 
effect  of  lamp  heating  in  the  ASM  Epsilon-1  reactor  for  use  in  high  fidelity  models  of  epitaxial  growth. 

Methodology 

The  analytical  approach  we  take  to  determine  the  heat  flux  spatial  profiles  is  based  on  the  concept  of 
view  factor  [37,  45]  which  describes  the  radiation  exchange  between  two  or  more  surfaces  separated  by  a 
nonparticipating  medium  that  does  not  absorb,  emit,  or  scatter  radiation.  The  view  factor  between  two 
surfaces  represents  the  fraction  of  radiative  energy  leaving  one  surface  that  strikes  the  other  surface  directly. 

In  this  method,  the  geometry  of  the  chamber,  including  location  and  shape  of  lamps,  susceptor,  reflectors, 
and  possibly  other  apparatus,  is  what  mainly  determines  the  form  of  the  resulting  flux  profiles.  This  geometric 
approach  was  adopted  in  [16],  where  the  authors  consider  only  a  two-dimensional  slice  of  the  chamber 
geometry,  and  includes  the  effect  of  reflectors  behind  the  lamp  banks.  There,  the  two-dimensional  approach 
was  reasonable,  perhaps,  since  the  lamp  arrangement  in  the  reactor  under  consideration  was  axisymmetric 
about  the  wafer  center.  This  situation  is,  however,  not  the  case  in  the  ASM  Epsilon-1  reactor.  Hence,  our 
analysis  is  similar  to  that  used  in  [10],  where  the  authors  consider  the  chamber  geometry  from  a  three- 
dimensional  point  of  view.  However,  in  that  paper,  as  in  this  paper,  the  effect  of  reflectors  is  not  included. 

Scope 

In  this  report,  radiant  flux  profiles  for  both  linear  and  spot  lamps  are  determined.  In  the  actual  reactor, 
the  internal  surface  of  the  chamber  lid  is  gold  plated  to  reflect  infrared  rays  from  the  linear  lamps,  and  the 
spot  lamps  are  placed  in  gold  plated  parabolic  reflectors.  However,  in  this  report,  we  do  not  consider  the 
effect  of  reflections  on  the  lamp  heating  of  the  wafer. 

Finally,  the  literature  indicates  that  “virtual  images” ,  or  radiation  from  the  heated  wafer  to  the  reflectors 
and  chamber  walls  which  is  reflected  back  to  the  wafer,  will  cause  additional  radiative  effects.  These  effects 
are  not  included  in  the  analysis  here. 

Assumptions 

We  shall  consider  all  surfaces  to  be  diffuse  reflectors  and  diffuse  emitters.  Radiant  intensity  from  lamps 
is  assumed  to  be  independent  of  direction  and  constant  across  the  length  of  the  lamp.  We  assume  that  the 
quartz  walls  and  the  process  gases  transmit  heat  radiation  from  the  lamps  perfectly  at  the  wavelengths  of 
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Radial  Position  Lamp 

On  Wafer  (inches)  Position 


-5.00 


1 1 .5  inches 

Figure  13:  Geometry  (top  view)  of  upper  lamp  array:  radial  position  of  each  lamp  is  given  in  inches  from 
center;  lamps  are  identified  by  five  uniquely  distinguishable  positions,  numbered  1  through  5. 


interest.  Furthermore,  we  assume  that  the  path  from  lamps  to  wafer  (or  lamps  to  susceptor)  is  completely 
free  of  any  other  obstacles. 

Relevant  Geometry 

Figure  13  shows  a  schematic  of  the  upper  lamp  array  superimposed  over  the  susceptor  and  wafer,  which 
is  based  on  a  description  and  diagram  provided  in  [5].  For  computational  purposes,  we  consider  each  linear 
lamp  to  be  a  straight  line  segment  of  length  11.5  inches  with  the  array  consisting  of  parallel  equally  spaced 
lamps.  The  array  begins  directly  above  the  susceptor  edge,  5.0  inches  (horizontally)  from  the  susceptor 
center.  The  distance  between  neighboring  parallel  lamps  in  the  array  is  1.25  inches.  The  vertical  distance 
between  wafer  and  upper  lamp  array  is  2.25  inches,  and  the  vertical  distance  between  wafer  and  lower  lamp 
array  is  3.50  inches.  We  note  that  the  distances  given  are  estimates  based  on  crude  measurements  taken  on 
the  reactor  itself. 

There  are  five  lamp  positions  for  the  linear  lamps  that  can  be  uniquely  differentiated  from  the  others. 
This  is  due  to  the  wafer  rotation.  For  example,  two  linear  lamps  equally  distant  from  the  center  linear  lamp 
have  an  identical  irradiating  effect  on  the  wafer  surface.  The  five  lamp  positions  are  numbered  1  through  5. 
The  spot  lamps  have  their  own  unique  geometry  and  are  analyzed  separately  later. 

The  source  of  radiation  for  each  lamp  is  a  tungsten  filament,  which  we  assume  to  be  a  straight  line 
segment  stretching  the  length  of  the  lamp.  Figure  14  shows  the  geometry  used  to  perform  the  analysis.  We 
assume  that  for  each  filament  the  radiant  intensity  is  independent  of  direction  and  constant  across  the  length 
of  the  filament. 

For  each  point  on  the  wafer,  w,  there  is  an  irradiance  contribution  from  each  point  on  the  filament,  /, 
depending  upon  the  distance  between  them,  d  =  \\w  —  /||,  the  angle  0W  formed  by  the  vector  w  —  f  and  the 
vector  nw  normal  to  the  wafer  surface,  and  the  vertical  distance,  h,  from  wafer  surface  to  filament.  Note  that 
on  the  filament  diagram  the  endpoint  values  are  /i  =  —5.75  inches  and  =  5.75  inches,  and  the  vertical 
distance  h  is  either  2.25  inches  for  the  upper  array  or  3.50  inches  for  the  lower  array. 

Heat  Flux  Calculation 

For  the  derivation  of  the  expression  for  heat  flux  radiant  power  per  unit  area  on  the  wafer  surface,  we 
adopt  the  notation  used  in  [37].  The  rate  of  radiative  energy  dQf  leaving  a  differential  surface  area  dAf 
(containing  the  point  /)  on  the  filament  that  strikes  a  differential  surface  area  dAw  (containing  the  point  w) 
on  the  wafer  surface  is  given  by 

dQf  =  dAf  If  cos{6f)  dcjfw  (29) 

where  If  is  the  intensity  of  radiative  energy  leaving  dAf  in  all  directions  in  hemispherical  space  (in  dimensions 
of  Watts  per  unit  area  per  steradian) ,  0 f  is  the  angle  formed  by  the  vector  w  —  f  and  a  vector  normal  to 
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Figure  14:  Geometry  for  view  factor  analysis  used  to  calculate  heat  flux  intensity 


profiles  for  linear  lamps. 


dAf ,  and  dujfw  is  the  solid  angle  subtended  by  dAw  from  /  given  by 


duJ  fw  — 


dAw  cos(Ow ) 

I2 


Substituting  (30)  into  (29)  yields 


dQf=dAfIf™ 


(30) 


(31) 


Now,  the  rate  of  radiation  energy  Qf  leaving  the  surface  element  dAf  on  the  filament  in  all  directions  over 
hemispherical  space  is  [37] 

Qf=itIfdAf.  (32) 

The  elemental  view  factor  dFdAf-dAw  is  defined  as  the  ratio  of  the  radiative  energy  leaving  dAf  that  strikes 
dAw  directly  to  the  radiative  energy  leaving  dAf  in  all  directions  into  the  hemispherical  space.  Thus,  we 
divide  (31)  by  (32)  to  give  the  view  factor 


dFdAf-dAw 


dQ  f  cos  ( 6 f )  cos  ( 0W )  dAw 

Qf  ird2 


(33) 


Since  we  are  assuming  that  filament  radiant  intensity  is  independent  of  direction,  we  take  Of  =  0  independent 
of  filament  position  /  so  that  cos  (Of)  =  1  and 


dFdAf-dAw 


cos(0w)  dAw 
nd2 


(34) 


We  are  interested  in  the  radiative  energy  illuminating  a  differential  area  on  the  wafer  due  to  the  entire 
filament.  To  compute  the  appropriate  view  factor,  Ff_dAw ,  we  average  (34)  across  the  length  of  the  filament 


Ff-dAw 


dAw  fh  cos(0w ) 
I/1-/2I4  ncP  h 


Finally,  we  observe  that 


for  the  given  geometry,  so  that 


cos(0w) 


h 

d 


Ff-dAw 


cL4w  fh  _h_ 

l/i  —  /2I  Jfl  nd3 


where  we  recall  that  d  =  \\w  —  /||. 


(35) 


(36) 
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To  determine  the  radiant  heat  flux  profile  for  a  given  lamp,  the  view  factor  Ff_dAw  must  be  computed 
for  each  differential  area  dAw  on  the  wafer  surface.  For  practical  purposes,  we  discretize  the  wafer  surface 
by  choosing  a  cylindrical  grid  of  wafer  points  w  =  (r,  (/)).  We  then  assume  that  the  differential  area,  dAw , 
is  constant  for  all  wafer  points  w.  Thus,  (36)  yields  the  view  factor  function  Ff_dAw  (4  <p)  which  gives  the 
fraction  of  radiative  energy  leaving  the  given  lamp  filament  that  strikes  the  given  wafer  point  w  =  (r,  q b) 
directly. 

Now,  we  let  Pf  denote  the  radiant  power  supplied  by  the  filament,  so  that  Pf  /dAw  gives  the  radiant  heat 
flux  intensity  striking  the  differential  area  dAw.  The  radiant  heat  flux  intensity  profile  of  the  illumination 
due  to  the  lamp  filament  is  then  given  by 


qf(r,4>)  =  Ff-dAw{r,4)) 

M 

7T  |/l  -  J 2|  Jfx 


dAu 

h 


f  -d_  df 

1-/2 1  4  ^M)3  1 


(37) 

(38) 


where  the  value  we  use  for  Pf  is  provided  by  the  manufacturer.  In  the  case  of  the  ASM  Epsilon-1  reactor, 
the  linear  lamps  supply  6000  Watts  and  the  spot  lamps  supply  1000  Watts. 

Since  the  wafer  rotates  at  a  uniform  rate,  this  function  is  averaged  over  the  circle  (i.e.,  0  <  <  2ir)  at 

each  radial  position  r  on  the  wafer  top  surface 


qf(r)  = 


dAv 


1  f27T 

2tt  J0 


Ffw(r ,  <t>)d<f> 


to  give  the  heat  flux  profile 


<lf(r) 


Pfh 


2tt2  |/i  -  /2| 


n/2 

.1 


d(r,  (f>)' 


■  df  d</). 


(39) 


(40) 


A  similar  analysis  was  performed  for  the  spot  lamps,  except  that  each  spot  lamp  was  considered  to  be  a 
point  source  of  radiant  energy,  thus  simplifying  the  analysis  significantly. 

The  computational  procedure  was  performed  for  each  of  the  five  different  linear  lamp  positions  for  both 
upper  and  lower  arrays,  and  for  the  spot  lamps.  Using  the  resulting  heat  flux  intensity  spatial  profiles, 
we  can  then  compute  the  desired  profiles  for  each  of  the  ten  lamp  groups,  and  the  four  heat  zones  of  the 
Epsilon-1  reactor  by  appropriately  combining  the  profiles  determined  from  the  individual  lamps. 

Results 

Here  we  discuss  some  results  of  the  analysis.  Note  that  in  what  follows,  the  term  “wafer  surface” 
may  represent  the  top  surface  of  the  wafer  and  exposed  susceptor,  or  the  bottom  surface  of  the  susceptor, 
depending  upon  the  lamp  group  being  considered. 

Individual  Lamps 

Figures  15,  16,  and  17  show  the  heat  flux  irradiated  on  the  wafer  surface  by  lamps  in  positions  1,  2,  3, 
and  4  of  the  upper  and  lower  array,  respectively,  and  the  spot  lamp  position.  As  expected,  points  on  the  wafer 
surface  directly  under  (or  over)  the  lamp  filament  receive  the  most  intense  illumination,  i.e.  the  maximum 
flux  value.  Intensities  are  greater  in  magnitude  for  lamps  in  the  upper  array  since  it  is  physically  closer  to 
the  wafer  than  the  lower  array  and  spot  lamps.  Spot  lamps  have  lower  flux  intensities  than  linear  lamps  due 
to  the  smaller  supplied  power. 

To  account  for  wafer  rotation,  the  flux  intensity  profiles  are  averaged  around  360  degrees  resulting  in 
profiles  that  are  a  function  of  radial  position  only.  Figures  18  and  19  show  the  heat  flux  profiles,  after 
averaging,  for  each  of  the  individual  lamp  positions.  Observe  that  as  expected  the  lamp  position  directly 
over  (or  under)  the  wafer  center  irradiates  the  wafer  center  with  greater  intensity  than  the  other  positions. 
Lamp  positions  closer  to  the  edge  irradiate  the  edge  with  greater  intensity  than  they  irradiate  the  center. 

Lamp  Groups 


Figure  20  shows  the  heat  flux  profiles  for  each  of  the  ten  lamp  groups. 
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Linear  Lamp  -  Upper  Array  -  Pos  1  Linear  Lamp  -  Upper  Array  -  Pos  2 


Figure  15:  Heat  flux  intensity  profiles  for  linear  lamps  in  upper  array:  flux  intensity  (W/cm2)  versus  position 
in  two  dimensions.  Upper  left:  Position  1;  Upper  right:  Position  2;  Lower  Left:  Position  3;  Lower  right: 
Position  4. 


Linear  Lamp  -  Lower  Array  -  Pos  1  Linear  Lamp  -  Lower  Array  -  Pos  2 


Figure  16:  Heat  flux  intensity  profiles  for  linear  lamps  in  lower  array:  flux  intensity  (W/cm2)  versus  position 
in  two  dimensions.  Upper  left:  Position  1;  Upper  right:  Position  2;  Lower  Left:  Position  3;  Lower  right: 
Position  4. 
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Spot  Lamp  Heat  Flux  Intensity 


Figure  17:  Heat  flux  intensity  profile  for  spot  lamp:  flux  intensity  (W/cm2)  versus  position  in  two  dimensions. 


Lower  Lamp  Array  -  Individual  Lamps  -  Heat  Flux  Intensity  Profiles 


Figure  18:  Heat  flux  intensity  profiles  for  individual  lamps  in  the  lower  array:  flux  intensity  (W/cm2)  versus 
radial  position  for  the  five  uniquely  distinguishable  linear  lamp  positions  and  the  spot  lamp  position. 
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Upper  Lamp  Array  -  Individual  Lamps  -  Heat  Flux  Intensity  Profiles 


Figure  19:  Heat  flux  intensity  profiles  for  individual  lamps  in  the  upper  array:  flux  intensity  (W/cm2)  versus 
radial  position  for  the  five  uniquely  distinguishable  linear  lamp  positions. 
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Figure  20:  Heat  flux  intensity  profiles  for  ASM  Epsilon-1  lamp  groups:  flux  intensity  (W/cm2)  versus  radial 
position  for  the  ten  lamp  groups. 
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Heat  Flux  Intensity  Profiles  for  ASM  Reactor  Lamp  Zones 
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Figure  21:  Heat  flux  intensity  profiles  for  ASM  Epsilon-1  heat  zones:  flux  intensity  (W/cm2)  versus  radial 
position  for  the  four  heat  zones. 


Heat  Zones 

Figure  21  shows  the  heat  flux  profiles  for  the  four  heating  zones  -  center,  front,  rear,  and  side.  The 
flux  intensity  for  the  center  zone  is  significantly  greater  than  for  the  others,  indicating  that  it  will  have  the 
greatest  heating  effect.  Observe  that  profiles  for  front  and  rear  zones  are  identical  due  to  the  symmetry 
assumptions  and  the  way  in  which  the  individual  lamps  are  organized  to  form  the  zones. 

Deficiencies  and  Improvements 

There  are  several  deficiencies  in  the  lamp  heating  model  as  it  currently  stands.  Improvements  depend 
on  gathering  of  more  accurate  chamber  geometry  data,  determination  of  lamp  and  reflector  properties,  and 
inclusion  of  more  complicated  effects. 

Reflectors 

The  internal  surface  of  the  chamber  lid  is  gold  plated  to  reflect  infrared  rays  from  the  lamps.  In  addition, 
the  spot  lamps  are  placed  in  gold  plated  parabolic  reflectors.  The  effects  of  reflections  must  be  included  in 
the  model.  In  order  to  compute  these  effects,  the  properties  and  geometry  of  the  reflectors  must  be  known. 

Accurate  Geometry 

The  analysis  has  proceeded  based  on  reactor  manual  diagrams  which  do  not  provide  exact  measurements 
for  the  chamber  geometry.  Accurate  dimensions  must  be  provided  to  improve  the  model. 

Lamp  Intensity 

The  lamp  heat  flux  intensity  profiles  have  been  determined  using  the  power  ratings  listed  for  the  lamps 
in  the  reactor  manual.  These  ratings  are  6  kW  for  the  linear  lamps  and  1  kW  for  the  spot  lamps.  However, 
we  do  not  know  the  efficiency  of  the  lamps  in  terms  of  proportion  of  power  consumed  to  power  supplied 
for  heating  the  wafer.  Thus,  the  profiles  may  be  off  by  a  scale  factor,  the  magnitude  of  which  cannot  be 
determined  exactly  without  more  information. 

Physical  Obstacles  and  Apparatus 
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Apparatus  or  other  physical  obstacles  in  the  reactor  may  have  an  effect  on  the  lamp  heating  by  blocking 
radiation  from  certain  directions.  Currently  these  are  not  considered.  If  such  obstacles  exist,  the  analysis 
will  be  more  complicated. 

Virtual  Images 

The  wafer  itself  will  irradiate  the  reflector  on  the  chamber  lid,  creating  virtual  images  of  the  wafer  which 
will  reflect  heat  flux  back  onto  the  wafer  surface.  These  effects  may  be  significant. 

2.3.5  Chemical  Reaction  Kinetics 

The  model  we  use  for  deposition  kinetics  is  based  on  the  Arrhenius  relationship  for  reaction  rate  described  in 
[20,  39,  42,  51].  We  assume  equilibrium  of  silane  and  hydrogen  adsorption  and  desorption,  and  that  the  rate 
of  silicon  deposition  is  limited  by  the  rate  at  which  surface  adsorbed  silane  decomposes  to  produce  silicon 
and  hydrogen  in  the  reaction 

SiH4  — »  Si  +  2H2.  (41) 

The  reaction  rate  at  each  position  on  the  wafer  surface  is  dependent  on  the  mole  fraction  of  silane,  XsiH4, 
and  the  wafer  temperature,  TWl  at  that  position.  It  is  given  by 

Rsi(Tw,XsiH4)  =  ks(Tw)X  siH4  (42) 

where  ks(Tw )  is  the  reaction  rate  parameter  given  by 

ks(Tw)  =  k0 exp  (43) 

and  ko  is  the  Arrhenius  coefficient,  Ea  is  the  activation  energy,  and  Rg  is  the  gas  constant. 

Equivalently,  we  can  replace  the  silane  mole  fraction  dependence  in  (42)  with  dependence  on  silane 
concentration,  CsiH4,  or  silane  partial  pressure,  PsiH4- 


Rsi(TW:  CsiH4)  =  k's(Tw)  Csm4  (44) 

Rsi(Tw,Psm4)  =  k"(Tw)Psm4  (45) 

where  in  each  case  the  value  and  units  of  the  Arrhenius  coefficient  ko  in  (43)  is  appropriately  converted  to  ac¬ 
count  for  the  change  of  units.  However,  for  our  purposes,  it  is  most  convenient  to  use  the  rate  expression  (42) 
since  the  mole  fraction  Xsih4  is  dimensionless  and  hence  all  units  can  be  accounted  for  in  ko- 

Values  for  the  activation  energy  Ea  and  the  Arrhenius  coefficient  ko  depend  on  the  process  operating 
conditions,  the  source  gases  used,  and  the  reactor  itself.  The  molecular  weight  and  mass  density  of  the 
deposited  material,  silicon  in  this  case,  is  also  incorporated  into  ko-  We  experimentally  determined  the 
values  for  ko  and  Ea  as  will  be  described  in  Section  2.4.1.  The  values  are  given  in  Appendix  A.  Note  that 
the  value  for  Ea  is  in  the  range  of  1.5-1. 7  eV  and  the  ratio  Ea/Rg  is  in  the  range  of  18000-20000,  both  of 
which  lie  within  the  range  of  values  typically  found  for  a  variety  of  operating  conditions  (see,  e.g.,  [39,  51]). 
Reaction  rate  as  a  function  of  time  and  radial  position  on  the  wafer  surface  is  given  by 

Rsi(t,r)  =  ko  exp  (fl^;r))  Xsm,(t,r) 

so  that 

af  (*,’’)  =  ■Rs‘(<,’')  =  toexp \RgTu,(t,r)  J  <47) 

describes  the  time  evolution  of  the  deposition  thickness  profile  h(t,  r).  Equation  (47)  is  numerically  inte¬ 
grated  in  conjunction  with  equation  (19),  the  evolution  equation  for  the  wafer  temperature  field  Tw(t,  r),  to 
determine  the  thickness  profile  of  deposited  silicon  on  the  wafer  surface.  The  mole  fraction  profile  Xsih4  ( t ,  r) 
can  be  determined  from  the  models  for  gas  flow  and  transport  of  chemical  species  in  the  reactor.  We  can 
also  simplify  the  situation  and  assume  a  pre-determined  constant  silane  mole  fraction  profile  Xsm4(r),  e.g., 
a  uniform  value  equal  to  the  average  silane  mole  fraction  found  during  a  particular  deposition  run. 


(46) 
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When  dimensionless  variables  are  used  to  avoid  scaling  problems  in  computational  work,  we  use  the 
conversions 


h  — )> 


t 


t 


href  T 

so  that  the  time  evolution  of  film  thickness  is  given  by 
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are  dimensionless  constants,  values  of  which  are  given  in  Appendix  A. 


2.3.6  Preliminary  Results 

Several  simulations  were  run  to  demonstrate  the  predictive  capabilities  of  the  wafer  heat  transfer  and  chemical 
kinetics  models.  These  simulations  are  tests  to  determine  if  the  models  yield  physically  reasonable  results.  In 
order  to  validate  time  evolution  of  the  wafer  temperature  field,  further  experimentation  will  be  necessary.  For 
example,  we  could  record  the  wafer  temperature  profile  as  a  function  of  time  using  an  array  of  thermocouples, 
and  compare  with  simulation  data.  Time  evolution  of  deposition  thickness  is  validated  by  comparing  with 
the  experimental  results  presented  in  Section  2.4.1. 

Wafer  Temperature  Ramp  Up 

In  this  simulation  we  start  the  wafer  at  a  uniform  room  temperature  and  apply  25%  power  to  all  lamp 
zones  for  one  minute.  Time  evolution  of  wafer  temperature  is  shown  in  Figure  22.  The  wafer  temperature 
increases  nommiformly  across  the  surface  due  to  the  greater  intensity  of  the  center  lamp  zone. 

During  ramp-up,  lamp  heating  dominates  the  heat  transfer  mechanisms  until  radiative  losses  become 
more  apparent  as  the  temperature  increases  beyond  ambient  chamber  wall  temperature.  Effects  of  conduction 
do  not  appear  since  the  temperature  profiles  are  relatively  uniform  throughout  and  the  intensity  of  radiative 
transfer  is  relatively  high.  The  average  temperature  levels  off  at  approximately  900  C. 

Wafer  Temperature  Ramp  Down 

In  this  simulation  we  start  the  wafer  at  a  uniform  processing  temperature  of  725  C  and  turn  off  all  lamp 
groups  for  5  minutes.  Time  evolution  of  wafer  temperature  is  shown  in  Figure  23.  The  wafer  temperature 
drifts  down  to  approximately  230  C  due  to  radiative  and  convective  losses  from  the  wafer.  Note  that  the 
temperature  profile  has  a  stable  equilibrium,  settling  to  a  spatially  uniform  profile  Tw(t,r )  — >  Tconst  for  all 
r  as  t  — >•  oo.  Here,  Tconst  is  equal  to  the  ambient  chamber  wall  temperature  Tc  when  convective  losses  are 
ignored,  and  lower  than  Tc  when  convective  losses  are  included.  Also,  note  that  the  temperature  changes 
more  rapidly  than  an  exponential  decrease  due  to  the  nonlinearity  of  the  radiative  effect.  Heat  diffusion  in 
the  solid  has  no  effect  in  this  simulation  since  the  temperature  remains  uniform  throughout. 

Temperature  Behavior  In  Response  To  Disturbances 

This  simulation  demonstrates  what  happens  when  the  temperature  profile  is  locally  perturbed  from 
a  uniform  profile.  Lamp  heating  is  turned  off  in  this  simulation  so  that  we  can  observe  the  effects  of  heat 
diffusion.  We  inject  a  temperature  spike  half-way  between  the  wafer  center  and  edge  and  record  temperature 
profiles  at  intervals  of  1.2  seconds  over  12  seconds.  The  results  are  shown  in  Figure  24.  Heat  conduction  in 
the  solid  wafer  causes  the  expected  smoothing  effect.  This  behavior  is  consistent  with  what  we  expect  to 
happen  under  our  modeling  assumptions,  i.e.,  the  wafer  thermal  dynamics  are  globally  asymptotically  stable 
about  a  uniform  spatial  profile  slightly  below  chamber  wall  temperature  (slightly  below  due  to  convective 
losses). 

Film  Growth  Under  Uniform  Temperature  Conditions 

We  demonstrate  film  growth  under  uniform  wafer  surface  temperature  conditions  by  simulating  the  film 
growth  under  uniform  temperature  profiles  of  700  C,  725  C,  750  C,  and  775  C  over  a  5  minute  deposition 
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Wafer  Temperature  -  All  Lamp  Zones  25%  Power 


Figure  22:  Time  evolution  of  temperature  profile  across  wafer  surface:  all  lamp  zones  at  25%  power  over  1 
minute. 


Wafer  Temperature  -  All  Lamp  Zones  Off 


Figure  23:  Time  evolution  of  temperature  profile  across  wafer  surface:  all  lamp  zones  at  0%  power  over  5 
minutes.  Note:  For  visual  clarity  axes  are  reversed  from  that  in  Figure  22. 
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Disturbance  Response  -  Transient  Decay  -  No  Heating 


Figure  24:  Response  to  temperature  disturbance  with  all  lamp  groups  at  0%  power  over  12  seconds. 


period.  Results  are  shown  in  Figure  25.  The  strong  temperature  dependence  of  the  growth  dynamics  are 
apparent.  These  simulation  results  match  the  experimental  results  from  growth  experiments  described  in 
Section  2.4.1. 

Film  Growth  Under  Nonuniform  Temperature  Conditions 

Here  we  simulate  growth  under  a  spatial  temperature  profile  that  linearly  increases  from  725.0  C  at  the 
center  to  725.7  C  at  the  edge.  This  0.7  degree  C  difference  amounts  to  a  1.2%  variation  in  wafer  temperature 
from  center  to  edge.  Growth  is  once  again  simulated  over  a  5  minute  deposition  period.  Results  are  shown  in 
Figure  26.  The  variation  in  temperature  is  amplified  into  a  19  Angstrom  or  26%  variation  in  thickness  from 
center  to  edge  at  the  end  of  the  5  minute  period.  This  clearly  demonstrates  the  importance  of  temperature 
control  for  reduced  deposition  nonuniformity. 

2.4  Experimental  Validation 

Andrew  Newman  of  ISR  and  Paul  Brabant  of  Northrop  Grumman  ESSD  conducted  thin  film  deposition  and 
lamp  heating  experiments  using  the  ASM  Epsilon-1  RTCVD  reactor  on  site  at  Northrop  Grumman  ESSD 
from  May  13,  1997  through  May  15,  1997.  Additional  measurements  were  taken  on  July  27,  1997. 

The  objectives  of  these  experiments  were  to 

•  validate  an  assumed  Arrhenius  relationship  between  wafer  temperature  and  deposition  rate,  and  de¬ 
termine  the  unknown  parameters  of  this  relationship  under  typical  operating  conditions; 

•  validate  the  correctness  of  analytically  determined  lamp  heating  models  described  earlier  in  this  report; 

•  measure  physical  dimensions  of  various  parts  of  the  reactor  in  order  to  determine  model  parameters; 
and 

•  familiarize  the  ISR  student  participant  with  the  operation  and  capabilities  of  the  reactor  and  various 
other  tools  used  at  the  Northrop  Grumman  ESSD  facility. 

The  ASM  Epsilon-1  reactor  at  Northrop  Grumman  ESSD  is  a  production  tool  in  regular  use.  The  authors 
are  grateful  for  the  time  and  resources  that  were  made  available  for  experimentation. 
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Temperature  =  700  C  Temperature  =  725  C 


Temperature  =  750  C  Temperature  =  775  C 


Figure  25:  Film  thickness  (Angstroms)  versus  time  (minutes)  for  a  range  of  wafer  temperatures  in  the 
thermally  activated  regime. 


Film  Thickness  -  Nonuniform  Temperature  Profile 


Figure  26:  Film  thickness  profiles  at  various  instants  of  time  under  nonuniform  wafer  temperature  conditions. 
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2.4.1  Deposition  Kinetics 

One  objective  of  our  experiments  was  to  validate  an  assumed  Arrhenius  relationship  between  wafer  tempera¬ 
ture  and  deposition  rate,  and  determine  the  unknown  parameters  of  this  relationship  under  typical  operating 
conditions.  Once  a  relationship  between  wafer  temperature  and  deposition  rate  was  established,  it  could  then 
be  used  in  the  model  for  deposition  kinetics.  Furthermore,  it  has  also  been  used  in  the  validation  analysis  for 
lamp  heating  models,  where  heat  flux  intensity  profiles  have  been  estimated  using  wafer  temperature  data 
which  in  turn  has  been  inferred  from  film  thickness  data.  This  analysis  is  described  later  in  Section  2.4.2. 

Although  the  overall  modeling  project  focuses  on  epitaxial  growth,  polycrystalline  silicon  was  deposited 
instead  of  epitaxial  silicon.  This  is  because  under  conditions  for  thermally  activated  growth,  deposition  rates 
are  relatively  slow,  making  impractical  the  growth  of  films  with  thicknesses  greater  than  one  micron.  For  epi¬ 
taxial  silicon,  available  instrumentation  at  Northrop  Grumman  (FTIR)  does  not  allow  for  the  measurement 
of  these  relatively  thin  (less  than  one  micron)  film  thicknesses.  On  the  other  hand,  polycrystalline  silicon 
film  thicknesses  between  5  Angstoms  and  1  micron  can  be  and  were  measured  by  the  available  ellipsometer 
and  nanospec.  Previous  work  by  Brabant  has  shown  growth  rates  for  epitaxial  and  polycrystalline  silicon  to 
be  approximately  the  same  [6] . 

Procedure 

Thin  films  of  polycrystalline  silicon  were  deposited  from  a  silane  precursor  over  a  five  minute  period  at 
reduced  pressure  (20  Torr).  Deposition  was  performed  under  a  combination  of  operating  conditions  consisting 
of  four  different  wafer  temperatures  and  three  different  silane  flow  rates  (and  hence  three  different  silane 
mole  fractions). 

Temperatures  were  set  in  the  surface-reaction  controlled  regime  so  that  growth  would  be  thermally 
activated.  This  regime  is  roughly  from  600  C  to  800  C  for  deposition  of  silicon  from  silane  gas.  We  chose 
the  following  wafer  temperatures  at  which  to  deposit  silicon:  650  C,  700  C,  725  C,  and  750  C. 

The  precursor  gas  mixture  was  silane  (SiFG)  diluted  at  2%  in  hydrogen  (H2).  Three  different  flow  rates 
were  used  for  this  2%  silane  in  hydrogen  precursor:  1.5  slm,  2.5  slm,  and  3.5  slm.  Considering  the  2% 
dilution,  these  three  flow  rates  correspond  to  30  seem,  50  seem,  and  70  seem  of  silane,  respectively.  The 
silane-hydrogen  precursor  was  again  diluted  in  20  slm  of  the  carrier  hydrogen  (H2)  gas.  Thus,  the  three  flow 
rates  correspond  to  three  mole  fractions  1.4  xlO-3,  2.2  xlO-3,  and  3.0  xlO-3,  respectively. 

The  reactor  was  operated  in  its  usual,  automatic  mode  (i.e.,  using  PID  control  loops  for  temperature 
regulation  and  wafer  rotation  for  uniformity),  using  pre-programmed  recipes.  Recipes  were  programmed  to 
set  chamber  pressure  at  20  Torr  and  to  deposit  silicon  from  silane  precursor  for  five  minutes  onto  the  bare 
silicon  wafers.  Film  thicknesses  were  measured  later  using  the  nanospec. 

Results 

We  attempted  twelve  deposition  experiments  -  one  for  each  combination  of  the  four  wafer  temperatures 
(650,  700,  725,  750  C)  and  three  silane  flow  rates  (30,  50,  70  seem).  Each  of  the  twelve  depositions  was 
performed  on  a  different  wafer.  At  650  C,  there  was  no  appreciable  deposition  for  any  of  the  flow  rates. 
Thus,  these  three  wafers  provided  no  data  for  analysis.  At  700  C  and  above,  enough  silicon  was  deposited 
so  that  measurements  could  be  taken.  Thickness  was  measured  using  the  nanospec  at  five  different  points 
on  the  wafer  surface.  These  five  points  are  shown  in  Figure  27.  In  the  case  where  temperature  was  700  C 
and  silane  flow  rate  was  30  seem,  film  thickness  was  less  than  100  Angstroms,  the  minimum  readable  by  the 
nanospec,  and  hence  is  recorded  as  less  than  100  Angstroms.  The  data  is  presented  in  Table  1. 

Deposition  kinetics  are  modeled  using  the  Arrhenius  relationship 

RSi  =  k0  exp  (^REr£^  ^SiH4  (48) 

where  Rsi  denotes  deposition  rate,  ko  denotes  the  pre-exponential  constant,  Ea  denotes  the  activation 
energy,  Rg  denotes  the  gas  constant,  Tw  denotes  the  wafer  temperature,  and  XsiH4  denotes  the  silane  mole 
fraction.  We  call  a  plot  of  the  logarithm  of  deposition  rate  versus  inverse  temperature  an  Arrhenius  plot.  The 
Arrhenius  plots  associated  with  the  data  we  collected  are  shown  in  Figure  28.  According  to  equation  (48), 
the  slope  of  the  Arrhenius  plot  gives  the  activation  energy  Ea  while  the  intercept  (along  with  knowledge  of 
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Silicon  Film  Thickness  and  Deposition  Rate 
As  Function  Of  Temperature  and  Silane  Flow  Rate 


Operating  Conditions 
Pressure  20  Torr 


Temperature  — 

700  C 

Silane  Flow  Rate  | 

30  seem 

50  seem 

70  seem 

Wafer  Position 

h(A) 

R  (A/min) 

h  (A) 

R  (A/min) 

h  (A) 

R  (A/min) 

1 

<  100 

<  20.0 

314 

62.8 

392 

78.4 

2 

<  100 

<  20.0 

337 

67.4 

393 

78.6 

3 

<  100 

<  20.0 

332 

66.4 

414 

82.8 

4 

<  100 

<  20.0 

332 

66.4 

391 

78.2 

5 

<  100 

<  20.0 

327 

65.4 

412 

82.4 

Avg 

<  100 

<  20.00 

328.4 

65.68 

400.4 

80.08 

Temperature  — 

725  C 

Silane  Flow  Rate  | 

30  seem 

50  seem 

70  seem 

Wafer  Position 

h  (A) 

R  (A/min) 

fa  (A) 

R  (A/min) 

fa  (A) 

R  (A/min) 

1 

365 

73.0 

526 

105.2 

684 

136.8 

2 

368 

73.6 

553 

110.6 

716 

143.2 

3 

365 

73.0 

539 

107.8 

693 

138.6 

4 

365 

73.0 

521 

104.2 

686 

137.2 

5 

365 

73.0 

526 

105.2 

689 

137.8 

Avg 

365.6 

73.12 

533.0 

106.60 

693.6 

138.72 

Temperature  — 

750  C 

Silane  Flow  Rate  | 

30  seem 

50  seem 

70  seem 

Wafer  Position 

fa  (A) 

R  (A/min) 

h  (A) 

R  (A/min) 

fa  (A) 

R  (A/min) 

1 

584 

116.8 

851 

170.2 

1072 

214.4 

2 

602 

120.4 

883 

176.6 

1116 

223.2 

3 

594 

118.8 

859 

171.8 

1082 

216.4 

4 

590 

118.0 

842 

168.4 

1070 

214.0 

5 

587 

117.4 

848 

169.6 

1077 

215.4 

Avg 

591.4 

118.28 

856.6 

171.32 

1083.4 

216.68 

Table  1:  Measured  silicon  film  thickness,  h  (Angstroms)  and  deposition  rate,  R  (Angstroms  per  minute): 
Five  minute  deposition;  three  wafer  temperatures;  three  silane  flow  rates. 
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Rear 


4  2  5  Side 

Center 


Front 


Figure  27:  Points  on  wafers  where  thickness  was  measured.  Rear  refers  to  side  of  wafer  closest  to  outlet 
(downstream);  front  closest  to  inlet  (upstream). 


the  silane  mole  fraction)  gives  the  pre-exponential  constant  fco-  Computed  parameters  are  given  in  Table  2. 
These  values  are  then  used  in  the  implementation  of  the  models  for  chemical  reaction  kinetics. 

The  activation  energies  range  from  1.57  eV  to  1.69  eV  depending  on  silane  mole  fraction.  This  range  is 
very  close  to  the  activation  energy  of  1.82  eV  determined  experimentally  by  the  manufacturer,  ASM,  Inc.  [29] 
In  addition,  the  pre-exponential  constants  range  from  3.8  x  108  to  1.85  x  109,  a  range  which  includes  the 
value  of  7.9  x  108  given  by  the  manufacturer. 

The  plots  of  deposition  rate  as  a  function  of  silane  flow  rate  and  silane  mole  fraction  are  shown  in 
Figures  29  and  30,  respectively.  As  expected,  the  relationship  is  nearly  linear,  with  slope  proportional  to 
the  exponential  of  inverse  temperature,  which  confirms  the  assumption  used  in  the  Arrhenius  model.  The 
relationship  between  deposition  rate  and  silane  mole  fraction  can  be  used  to  determine  the  effect  of  reactant 
depletion  across  the  wafer  surface  during  deposition. 

2.4.2  Lamp  Heating 

The  objective  of  lamp  heating  experiments  was  to  provide  a  basis  for  comparison  to  validate  the  analytically 
determined  heat  flux  intensity  profiles  described  in  Section  2.3.4.  In  particular,  the  goal  here  was  to  isolate 
individual  lamp  groups  in  order  to  determine  their  heating  effect  on  the  wafer  surface.  Recall  that  lamp 
heating  enters  the  wafer  heat  transfer  ODE  (19)  via  the  lamp  influence  matrix  B  in  the  control  input  term, 
which  was  determined  using  analytical  methods.  It  can  be  compared  with  experimentally  determined  lamp 
influence  functions  for  purposes  of  validation. 

Recall  that  for  each  individual  lamp  we  computed  a  spatial  profile  that  gives  a  value  of  heat  flux  intensity 
(in  W / cm2)  irradiating  each  point  on  the  wafer  surface,  where  points  are  defined  by  their  radial  and  azimuthal 
coordinates.  Then,  spatial  profiles  were  averaged  around  all  azimuthal  positions  to  simulate  wafer  rotation 
and  provide  1-dimensional  profiles  that  are  functions  of  radial  position  only.  Finally,  profiles  for  individual 
lamps  were  combined  appropriately  to  give  heat  flux  intensity  profiles  for  the  ten  lamp  groups  and  four  lamp 
zones  of  the  ASM  Epsilon-1  reactor. 

Analytically  determined  influence  functions  can  be  computed  using  an  arbitrarily  fine  spatial  resolution. 
The  only  cost  is  computing  time,  and  the  computations  require  relatively  little  time  on  a  high  performance 
workstation.  For  example,  computation  of  heat  flux  intensity  for  one  individual  lamp  on  a  cylindrical  grid 
with  3200  points  (100  radial,  32  azimuthal)  required  less  than  3  minutes.  As  we  shall  see,  the  procedure  that 
we  used  for  experimentally  determination  requires  a  physical  measurement  of  film  thickness  at  each  spatial 
point.  It  is  apparent  that  doing  this  for  a  grid  with  a  large  number  of  points  (say  more  than  100)  becomes 
impractical.  Therefore,  the  experimentally  determined  profiles  are  used  only  for  validation. 

Methodology 

The  following  methodology  was  used  to  empirically  determine  the  lamp  influence  functions.  Details  of 
the  experimental  procedure  are  described  later. 

1.  Polysilicon  is  deposited  on  a  non-rotating  wafer  for  a  fixed  period  of  time  r  with  lamp  group  i  manually 
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Arrhenius  Plot  for  Poly-Si  Deposition  Rate 


Figure  28:  Arrhenius  plots  for  silicon  deposition  from  silane  gas:  each  plot  represents  log  of  deposition  rate 
(microns  per  minute)  versus  inverse  absolute  temperature  for  one  of  the  three  silane  flow  rates  used. 


Parameters  For  Arrhenius  Relationship  Describing  Silicon  Deposition  Kinetics 


Rsi 


ko  exp 


~Ea  \ 

7?  T 


XsiH4 


Symbol 

Description 

Data 

Fmix 

Silane/Hydrogen  Mixture  Flow  Rate  (slm) 

1.5 

2.5 

3.5 

VsiH4 

Silane  Flow  Rate  (seem) 

30 

50 

70 

A^SiH4 

Silane  Mole  Fraction  (xlO-3) 

1.4 

2.2 

3.0 

Ea 

Activation  Energy  (eV) 

1.69 

1.67 

1.57 

Activation  Energy  (J/mol)  (xlO5) 

1.63 

1.61 

1.51 

Ea/  Rg 

Ratio  (K)  (xlO4) 

1.96 

1.94 

1.82 

k0 

Pre-exponential  Constant  (um/min)  (xlO9) 

1.85 

1.30 

0.38 

ko 

Pre-exponential  Constant  (cm/sec)  (xlO3) 

3.08 

2.16 

6.54 

Table  2:  Parameters  for  Arrhenius  relationship  describing  silicon  deposition  kinetics:  activation  energy  (eV) 
and  pre-exponential  constant  (microns  per  minute)  for  each  silane  flow  rate  used. 
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Deposition  Rate  As  Function  Of  Silane  Flow  Rate 


Figure  29:  Growth  rate  as  function  of  silane  flow  rate  for  each  of  three  temperatures  used. 


Deposition  Rate  As  Function  Of  Silane  Mole  Fraction 


Figure  30:  Growth  rate  as  function  of  silane  mole  fraction  for  each  of  three  temperatures  used. 
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set  to  power  setting  P.  Measuring  thickness  yields  a  thickness  function  /iq;p;T)(r,  9).  Growth  rate  is 
computed  as  R(i:p)(r,6)  =  h/r.  Averaging  azimuthally  gives  growth  rate  in  terms  of  radial  position 
%p)M- 

2.  The  Arrhenius  law  (46)  is  inverted  to  determine  temperature  as  a  function  of  radial  position 

%P)(r)  =  (Ea/Rg)  (Info  VsinJ  -  ]n(%P)(r)))_1  (49) 

3.  The  temperature  field  7q?p)(r)  is  substituted  into  the  evolution  equation  for  temperature  (19)  in  the 
wafer  heat  transfer  model.  By  applying  the  steady-state  condition  T  =  0  we  can  solve 

0  =  Ac  T  +  Ar  T4  +  Av  T  +  T  +  Bi  P  (50) 

for  the  discretized  influence  function  Bi(r). 

Experimental  Procedure  and  Results 

Isolation  of  the  individual  lamp  groups  was  achieved  by  operating  the  reactor  in  manual  mode,  i.e.,  with 
the  automatic  control  loops  for  temperature  regulation  turned  off.  In  manual  mode,  the  lamp  groups  are 
no  longer  organized  into  four  zones  for  the  purpose  of  temperature  control.  Instead,  each  of  the  ten  groups 
can  be  toggled  on  and  off  individually,  and  the  power  setting  of  each  (between  0%  and  100%)  can  be  set 
manually.  To  isolate  a  particular  lamp  group,  all  others  were  turned  off,  while  the  power  setting  for  the  lamp 
group  being  tested  is  set  manually  to  an  appropriate  level. 

The  wafer  was  heated  with  an  individual  lamp  group,  whose  power  setting  was  adjusted  manually,  until 
at  least  one  of  the  thermocouple  readings  reached  the  range  where  deposition  would  occur,  approximately 
700  C.  The  exact  temperature  readings  were  not  important  because  in  the  next  step  temperature  would  be 
inferred  from  thickness  data.  Then,  flow  of  silane  in  hydrogen  carrier  was  started.  Silicon  was  deposited 
for  five  minutes.  Wafer  rotation  was  turned  off  so  that  effects  of  asymmetry  would  appear  in  the  resulting 
deposition. 

This  procedure  was  followed  to  test  four  of  the  lamp  groups:  1,  8,  9,  and  10.  Lamp  group  1  is  in  the 
upper  lamp  array  and  radiates  directly  toward  the  top  center  of  the  wafer.  Using  lamp  group  1  alone,  we 
were  able  to  heat  the  wafer  to  a  temperature  sufficiently  high  for  deposition  to  occur  and  to  record  sufficient 
data  for  analysis.  Lamp  groups  8,  9,  and  10  are  in  the  lower  lamp  array  and  radiate  toward  the  bottom 
of  the  susceptor.  Due  to  conduction  and  losses  throughout  the  susceptor,  it  was  more  difficult  to  heat  the 
wafer  using  each  of  these  lamp  groups  alone.  Of  the  lamp  groups  we  isolated  in  the  lower  array,  only  lamp 
group  8  provided  enough  radiant  energy  to  heat  the  wafer  to  a  temperature  sufficiently  high  for  deposition 
to  occur.  However,  wafer  temperature  oscillated  and  was  highly  nonnniform  in  this  case,  causing  the  data 
to  be  unreliable.  We  focus  now  on  the  experiment  that  tested  lamp  group  1  from  which  reliable  data  was 
obtained. 

Lamp  group  1  was  isolated  and  set  to  45%  of  full  power  which  brought  the  center  thermocouple  reading 
to  740  C,  sufficiently  high  for  deposition  to  occur.  Silane  flow  rate  was  set  at  30  seem.  After  a  five  minute 
deposition  period,  the  wafer  was  removed  and  thickness  measurements  were  taken  at  100  points  on  the  wafer 
surface  as  shown  in  Figure  31. 

Figure  32  shows  two  views  of  the  resulting  polysilicon  film  thickness  profile.  Thermally  activated  growth 
using  the  isolated  lamp  group  1  produces  a  “hill”  of  polysilicon.  The  deposition  pattern  reaches  a  maximum 
in  a  line  across  the  wafer  center  parallel  to  the  lamps  in  lamp  group  1,  and  decreases  toward  the  wafer  edges. 
Qualitatively,  this  result  is  what  we  would  expect  given  the  geometry  of  this  particular  experimental  setup. 

The  thickness  data  is  then  used  to  compute  an  empirically  determined  heat  flux  intensity  profile  for  lamp 
group  1  as  outlined  previously.  Figure  33  shows  the  result  along  with  the  analytically  determined  profile 
for  purposes  of  comparison.  The  result  indicates  a  reasonable  agreement  between  the  analytical  model  and 
experimental  data. 

2.5  Feature  Scale  Models 

In  order  to  study  deposition  on  patterned  wafers,  the  macroscopic  level  equipment  and  process  model  will 
be  integrated  with  a  feature  scale  model.  This  model  should  describe  epitaxial  growth  of  thin  films  at  the 


36 


Front 


Figure  31:  Points  on  wafer  surface  where  thickness  measurements  were  taken  for  testing  of  lamp  group  1. 

microscopic  level,  i.e.,  based  on  the  physics  of  atom  motion  on  surfaces.  In  this  way,  we  can  study  how  the 
oxide  pattern  on  a  wafer  affects  the  depostion  process,  and  in  particular  how  pattern  pitch  affects  uniformity. 

The  motion  of  atoms  “raining  down”  on  the  surface  of  a  wafer  during  the  deposition  process  largely 
determines  the  arrangement  of  these  atoms  in  terms  of  crystal  structure  and  ultimately  the  quality  of  the 
film  for  its  intended  purpose  [27].  One  framework  for  studying  atomic  mechanisms  of  crystal  growth  is  the 
terrace-ledge-kink  model  introduced  by  Burton,  Cabrera,  and  Frank  [7].  The  general  idea  is  illustrated  in 
Figure  34.  In  this  framework,  arriving  atoms  land  on  the  surface  which  contains  terraces,  steps,  vacancies, 
and  kinks.  Atoms  can  evaporate,  move  on  the  terrace,  cross  over  steps,  nucleate  new  islands  on  any  terrace, 
or  become  incorporated  into  steps.  Transport  of  atoms  over  and  along  steps  is  the  controlling  factor  in 
epitaxial  growth  technology,  because  of  its  influence  on  the  film  morphology,  i.e.,  roughness  [27]. 

One  useful  model  for  our  purposes  here  is  the  “Eden  model”  for  growth  of  interfaces  relaxing  by  surface 
diffusion.  In  [13],  an  evolution  equation  for  film  growth  is  given  under  the  assumption  that  the  growth  rate 
depends  only  on  the  local  surface  morphology  through  rotational  invariants  (such  as  curvature), 

^  =  -K(V2)2h  +  uV2h  +  b(V/i)2  +  rj  (51) 

C/L  A 

where  h  =  h(r,  t )  is  the  height  (thickness)  profile,  r]  =  77 (r,  t )  is  the  statistical  noise  of  incoming  particle  flux, 
and  K  and  v  are  identified  from  physical  considerations. 

We  are  currently  investigating  the  established  body  of  work  in  this  area  in  the  hope  of  developing  an 
adequate  microscopic  model  for  purposes  of  modeling  epitaxial  growth  of  silicon  and  Si-Ge  on  wafers  with 
oxide  patterns. 

3  Model  Reduction 

Dynamical  system  models  have  been  presented  in  this  report  that  describe  the  time  evolution  of  various 
physical  phenomena  involved  in  RTCVD.  Model  reduction  deals  with  methods  for  reducing  the  dimensionality 
of  dynamical  system  models.  The  motivation  is  that  models  of  lower  dimension  are  less  complex  and  easier 
to  work  with  for  various  purposes  such  as  simulation,  optimization,  and  control.  In  this  section,  we  describe 
our  recent  work  in  model  reduction,  mainly  as  applied  to  the  models  for  wafer  heat  transfer  described  in  this 
report.  We  also  briefly  discuss  ongoing  work  to  extend  these  ideas  to  the  overall  reactor  model  including 
evolution  of  gas  temperature  and  chemical  species  transport. 
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Poly-Si  Thickness:  5  Min  Deposition  Using  Lamp  Group  1  At  45% 
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Poly-Si  Thickness:  5  Min  Deposition  Using  Lamp  Group  1  At  45% 


Figure  32:  Two  views  of  polysilicon  film  thickness  profile  resulting  from  5  minute  deposition  using  lamp 
group  1  at  45%  power  and  silane  flow  rate  of  30  seem.  Top  figure  shows  contour  map  where  colors/shades 
represent  thicknesses.  Bottom  figure  shows  3-dimensional  view  (“hill”). 
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Lamp  Group  #1  Heat  Flux  Intensity  Profiles 


Figure  33:  Experimentally  determined  heat  flux  intensity  profile  for  lamp  group  1  along  with  analytically 
determined  profile  for  purposes  of  comparison. 


Figure  34:  Terrace-ledge-kink  framework  for  describing  atom  motion  on  surfaces. 
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3.1  Motivation  and  Overview 


As  described  in  Section  2  of  this  report,  process  and  equipment  models  for  RTCVD  consist  mainly  of 
balance  equations  for  conservation  of  energy,  momentum,  and  mass,  along  with  equations  that  describe 
the  relevant  chemical  mechanisms.  In  their  continuum  form,  the  balance  equations  for  RTCVD  yield  a 
system  of  nonlinear  coupled  PDEs  with  associated  BCs  and  ICs.  Lumped  versions  of  the  equations  can 
be  obtained  using  an  appropriate  discretization  method,  e.g.,  finite-elements,  to  yield  a  system  of  coupled 
nonlinear  ODEs.  The  system  of  ODEs  can  be  decoupled  by  invoking  certain  simplifying  assumptions.  Even 
the  resulting  simplified  nonlinear  system  is  typically  of  relatively  high-order,  so  that  not  only  is  the  model 
computationally  demanding  for  simulation,  but  moreover  it  is  computationally  prohibitive  for  real-time 
control.  Thus,  the  motivation  for  reducing  the  model  order  is  apparent. 

A  general  approach  to  model  reduction  is  to  find  a  coordinate  transformation  of  the  state  space  under 
which  the  state  components  can  be  ranked  in  a  meaningful  way  in  terms  of  their  influence  on  the  system 
behavior.  Then,  state  components  of  the  transformed  system  with  relatively  small  influence  can  be  truncated 
without  substantially  degrading  the  correctness,  i.e.,  predictive  capability,  of  the  model.  We  note  that  for 
systems  evolving  on  lRn,  each  coordinate  transformation  can  be  identified  with  a  corresponding  set  of  basis 
n-vectors. 

To  illustrate  ideas  in  model  reduction,  we  focus  on  applying  model  reduction  methods  to  the  dynam¬ 
ical  system  ODE  model  describing  wafer  heat  transfer,  given  by  (19).  We  examine  two  model  reduction 
approaches  in  this  regard:  the  proper  orthogonal  decomposition  (POD)  and  the  method  of  balancing.  A 
comparison  of  the  effectiveness  of  the  two  approaches  is  presented  via  numerical  studies  using  the  wafer  heat 
transfer  model. 


3.2  Simplified  Wafer  Heat  Transfer  Model 


Recall  that  in  Section  2.3.3,  the  evolution  of  the  temperature  field  on  the  wafer  surface  was  given  by  the 
ODE 


Tw  —  Ac  Tw  +  Ar  T4  +  Av  Tw  +  T  +  5P 


(52) 


with  the  initial  condition 

TW(0)  =  TWO  (53) 

where  Tw0  represents  the  discretization  of  a  typical  inital  temperature  profile,  e.g.,  uniform  ambient  (ambient 
temperature  across  entire  surface). 

To  model  the  measurement  of  temperature  at  discrete  points  on  the  wafer  surface  via  thermocouples,  we 
augment  the  nonlinear  state  equation  (52)  with  the  linear  output  equation 


Ttc  =  CTw  (54) 

where  Ttc  is  a  p-vector  of  thermocouple  measurements,  and  C  is  a  p  x  n  matrix  with  entries  corresponding 
to  thermocouple  locations. 

As  described  in  Section  2.3  of  this  report,  the  various  parameters  in  (52)  and  (54)  are  derived  from  a 
detailed  analysis  of  the  ASM  Epsilon-1  reactor.  This  includes  the  lamp  heat  flux  intensity  profiles,  thermal 
parameters,  and  optical  parameters  used  in  the  model.  Under  our  modeling  assumptions,  the  reactor  has 
three  independent  lamp  actuators,  called  lamp  zones,  which  provide  m  =  3  independent  influence  functions. 
During  a  deposition  process,  there  are  no  thermocouples  on  the  wafer  surface.  However,  for  purposes  of  this 
study,  we  assume  that  there  are  thermocouples  placed  at  p  =  3  locations  on  the  wafer  surface:  center,  edge, 
and  midpoint  between  center  and  edge. 

Later  in  this  report  we  use  a  linearized  version  of  (52).  To  linearize,  first  observe  that 


x  —  Ac  x  T  Ar  (x  +  T)4  —  Ar  T4  T  Av  x  T  B  u  (55) 

has  an  equilibrium  point  at  x  =  0  and  is  equivalent  to  (52)  using  the  changes  of  variable  x  =  Tw  —  T  and 
u  =  P  —  Pss,  where  Pss  is  the  control  input  that  results  in  a  steady  state  temperature  field  of  Tw  =  T. 
Linearizing  (55)  about  the  origin  gives 

x  =  Ax  +  Bu  (56) 
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with 


A  —  Ac  +  Av  -f-  4  F 


(57) 


where 

W\ij  =  [Ar\ij  Tj 

and  x  and  u  are  translations  of  Tw  and  P,  respectively. 

As  before,  all  evolution  equations  are  numerically  integrated  using  a  fourth  and  fifth  order  Runga-Kutta 
algorithm,  and  the  number  of  states,  i.e.,  discretization  points,  for  the  original  model  and  the  linearized 
model  is  set  at  n  =  101. 

3.3  Proper  Orthogonal  Decomposition 

One  approach  to  finding  a  basis  for  the  desired  coordinate  transformation  is  application  of  the  POD,  which 
is  based  on  the  classical  Karhunen-Loeve  decomposition  of  a  stochastic  process  (see,  e.g.,  [14,  18,  35,  47]). 
The  POD  is  also  known  as  the  method  of  empirical  eigenfunctions,  or  principal  components  analysis  (PCA). 
The  POD  is  a  statistical  pattern  analysis  technique  for  finding  the  dominant  structures  in  an  ensemble  of 
spatially  distributed  data.  These  structures  can  be  used  as  an  orthogonal  basis  for  efficient  representation 
of  the  ensemble.  In  particular,  the  POD  basis  elements  are  the  orthogonal  eigenfunctions  of  the  two-point 
spatial  covariance  of  the  data  ensemble.  When  the  data  ensemble  consists  of  vectors  in  lRn,  the  POD  basis 
vectors  are  just  the  the  columns  of  the  matrix  V  in  the  singular  value  decomposition 

X  =  UEVT  (58) 

where  X  is  a  matrix  whose  columns  are  the  members  of  the  data  ensemble. 

For  the  case  of  a  dynamical  system  describing  a  flow,  e.g.,  a  heat  flow,  the  data  ensemble  is  typically 
time  series  data,  i.e.,  “snapshots”  of  the  flow  captured  at  equally  spaced  intervals  in  time.  The  time  series 
data  is  typically  generated  by  simulating  the  original  evolution  equations,  driven  by  an  ensemble  of  repre¬ 
sentative  control  inputs,  and  using  an  ensemble  of  representative  ICs.  The  POD  coordinate  transformation 
diagonalizes,  or  decouples,  the  covariance  of  the  time  series  data.  The  basis  elements  of  the  coordinate 
transformation  are  the  principal  axes  of  the  flow  which  generated  the  time  series  empirical  data.  Each  has 
a  corresponding  eigenvalue  (given  by  the  entries  on  the  diagonal  of  E  in  (58)),  which  provides  a  measure 
of  the  relative  “energy”,  i.e.,  mean  square  fluctuation,  associated  with  that  particular  direction  in  the  state 
space.  This  measure  can  also  be  interpreted  as  the  relative  amount  of  time  that  the  flow  spends  along  the 
corresponding  principal  axis.  Thus,  it  serves  as  a  well-defined  measure  of  the  influence  of  a  state  component 
on  the  system  behavior. 

The  POD  basis  is  optimal  from  the  points  of  view  of  data  compression  and  error  minimization.  Specifi¬ 
cally,  a  truncated  series  representation  of  the  data  has  a  smaller  mean  square  error  than  a  representation  by 
any  other  basis  of  the  same  dimension  (see,  e.g.,  [19]). 

Application 

The  attractive  properties  of  the  POD  have  led  to  success  in  applying  it  to  areas  such  as  turbulence 
modeling  (e.g.,  [19])  and  pattern  recognition  (e.g.,  [46]).  Recently,  much  work  has  been  done  to  study  its  use 
for  RTCVD  model  reduction  (e.g.,  [1,  2,  3,  4,  36,  50]).  We  have  performed  a  similar  study  toward  finding  a 
low-order  approximation  to  (52). 

To  generate  empirical  time  series  data,  i.e.,  snapshots  of  the  wafer  temperature  field,  the  system  (52-54) 
was  simulated  using  two  different  types  of  control  input  recipes.  (A  recipe  refers  to  a  function  giving  the 
lamp  power  setting  for  each  of  the  three  lamp  zones  at  each  instance  of  time.)  They  are  referred  to  as 
Ramp-Soak-Cool  (RSC)  and  Perturbation-Of-Constant  (POC). 

Control  Input  Recipe  Type  1  -  Ramp -Soak- Cool 

The  RSC  recipe  mimics  a  typical  processing  recipe  in  which  a  lamp  zone  power  setting  is  gradually  ramped 
up  from  zero  to  full  power,  maintained  at  full  power  for  a  specified  period  of  time,  and  then  gradually  ramped 
down  from  full  to  zero  power,  as  shown  in  Figure  35.  This  recipe  is  applied  to  one  of  the  lamp  zones,  while 
the  other  two  zones  are  held  at  zero  power.  The  simulation  is  then  repeated  using  RSC  individually  for  each 
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Ramp-Soak-Cool  Lamp  Power  Recipe 


Figure  35:  Lamp  power  settings  for  RSC  recipe. 


Radial  Position  (inches) 

Figure  36:  Snapshots  of  wafer  temperature  field  with  RSC  input  and  uniform  IC. 
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Figure  37:  Basis  elements  computed  using  POD  from  RSC  empirical  data. 


Lamp  Power  Settings  for  POC  Recipe 


Note:  Settings  are  constant  for  all  time. 


Recipe 

Zone  1 

Zone  2 

Zone  3 

Punif 

0.0798 

0.4265 

0.1965 

p{  1) 

0.0718 

0.4265 

0.1965 

p(2)  ' 

0.0878 

0.4265 

0.1965 

p(  3) 

0.0798 

0.3838 

0.1965 

p(4)  ‘ 

0.0798 

0.4691 

0.1965 

p(5) 

0.0798 

0.4265 

0.1768 

p(6)  ' 

0.0798 

0.4265 

0.2161 

Table  3:  Lamp  power  settings  for  POC  recipe. 


of  the  other  two  lamp  zones.  In  this  manner,  the  system  response  to  excitation  from  an  RSC  recipe  for  each 
of  the  three  lamp  zones  will  appear  in  the  time  series  data.  The  time  series  data  is  shown  in  Figure  36. 

The  entire  ensemble  (three  sets)  of  time  series  data  is  combined  and  arranged  into  a  data  matrix,  each 
column  of  which  represents  one  “snapshot”  of  the  wafer  temperature  field.  The  POD  basis  elements  and 
associated  eigenvalues,  or  relative  energy  values,  are  then  computed  via  SVD  and  ranked  according  to 
magnitude  of  relative  energy.  The  basis  elements  with  the  four  largest  eigenvalues  are  shown  in  Figure  37. 
Corresponding  relative  energy  values  are  contained  in  Table  4. 

Control  Input  Recipe  Type  2  -  Perturbation  Of  Constant 

The  POC  recipe  applies  small  perturbations  of  a  predetermined  set  of  constant  power  settings  which,  if 
left  unperturbed,  would  result  in  a  uniform  steady  state  temperature  field  of  1000K.  The  perturbations  are 
achieved  by  adjusting  the  power  setting  of  each  lamp  zone,  one  at  a  time,  first  to  110%  and  then  to  90%,  of 
the  original  setting.  This  results  in  6  different  control  recipes,  as  shown  in  Table  3.  Note  that  if  the  nominal 
constant  power  settings  are  used,  then  the  wafer  temperature  field  will  evolve  as  a  uniform  field  for  all  time. 
Thus,  the  perturbations  are  used  to  elicit  a  response  that  would  be  characterstic  of  the  system  behavior  in 
response  to  certain  types  of  disturbances. 

The  system  response  to  excitation  from  each  of  the  six  POC  recipes  is  sampled  and  combined  as  the  time 
series  data  for  computing  POD  basis  elements.  Time  series  snapshots  are  shown  in  Figure  38.  Once  again, 
POD  basis  elements  are  computed  and  ranked  by  corresponding  relative  energy  value.  The  basis  elements 
with  the  four  largest  eigenvalues  are  shown  in  Figure  37.  Corresponding  relative  energy  values  are  contained 
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Wafer  Temperature  Field  Snapshots  for  POC  Inputs 
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Figure  38:  Snapshots  of  wafer  temperature  field  with  POC  input  and  uniform  IC. 


Figure  39:  Basis  elements  computed  using  POD  from  POC  empirical  data. 


in  Table  4. 

Shortcomings 

It  is  clear  that  the  efficiency  of  a  basis  determined  via  the  POD  method  depends  strongly  on  how 
well  the  data  ensemble  captures  the  relevant  system  behavior.  This  leads  to  serious  shortcomings  for  model 
reduction  of  control  systems  with  inputs  and  outputs.  The  POD  basis  elements  will  be  sensitive  to  the  choice 
of  control  inputs,  ICs,  and  BCs  used  to  generate  the  empirical  time-series  data.  For  nonlinear  systems,  small 
perturbations  in  these  parameters  may  produce  qualitatively  different  system  responses.  In  addition,  the 
data  may  fail  to  capture  dynamical  effects  occuring  at  widely  differing  time  scales.  Usually,  these  issues 
are  ignored,  because  the  model  is  only  being  used  for  a  particular  purpose,  e.g.,  to  simulate  tracking  of 
one  particular  reference  trajectory.  In  that  case,  a  representative  set  of  control  inputs,  BCs,  and  ICs  is 
selected  for  generating  the  data  ensemble.  But  the  optimality  property  of  the  POD  basis,  and  the  predicitive 
capability  of  the  resulting  reduced  order  model,  may  be  localized  to  a  relatively  small  region  of  the  space 
of  allowable  inputs,  BCs,  and  ICs.  In  addition,  the  POD  approach  fails  to  consider  the  influence  of  state 
components  on  the  system  measurements,  or  outputs.  It  would  be  desirable  to  truncate  state  components 
that  have  small  influence  on  the  outputs,  i.e.  that  do  not  appear  in  our  measurements.  The  POD  method 
does  not  do  this,  thus  its  efficiency  may  be  diminished  when  constructing  “black-box”  input-output  models 
of  the  system. 
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3.4  Balancing 

In  response  to  concerns  regarding  the  POD  method,  we  considered  an  alternative  approach  based  on  the 
method  of  balanced  realizations  (see,  e.g.,  [12,  33]).  In  this  method,  a  coordinate  transformation  is  computed 
which  allows  state  components  to  be  ranked  according  to  their  influence  on  the  input-output  behavior  of 
the  system  as  measured  by  the  Hankel  norm  of  the  system,  i.e.,  the  gain  from  past  inputs  to  future  outputs. 
The  resulting  basis  for  the  linear  transformation  makes  the  transformed  realization  “equally  controllable 
and  observable,”  i.e.  balanced,  and  depends  only  upon  intrinsic  properties  of  the  original  model,  specifically 
controllability  and  observability,  embodied  in  its  evolution  equation  and  output  equation.  In  the  linear 
case,  explicit  bounds  can  be  computed  for  the  error  between  the  original  and  reduced  order  models.  These 
error  bounds  are  independent  of  any  particular  set  of  control  input  signals,  BCs,  or  ICs.  Although  explicit 
error  bounds  have  not  yet  been  found  for  the  nonlinear  case,  we  still  wish  to  exploit  the  property  that  the 
correctness  of  an  approximation  using  a  truncated  balanced  realization  does  not  depend  upon  generating 
a  representative  data  ensemble.  Furthermore,  the  balancing  method  emphasizes  state  components  that 
are  both  strongly  controllable  and  strongly  observable,  so  that  state  components  which  are  least  likely  to 
influence  the  measurements  are  truncated. 

The  theory  of  balancing  for  linear  systems  is  well  established.  See  Appendix  D  for  an  exposition  of  the 
main  ideas  and  results.  Computation  of  a  basis  for  the  balanced  realization  is  readily  accomplished  by  use  of 
matrix  operations  and  decompositions.  A  general  theory  and  procedure  for  balancing  of  a  class  of  nonlinear 
systems  has  been  presented  [43],  but  in  contrast  to  the  linear  case,  algorithms  for  performing  the  required 
computations  do  not  currently  exist.  We  note  that  development  of  such  algorithms  is  a  topic  of  current 
research  being  undertaken  by  the  authors.  However,  in  the  meantime,  application  of  the  balancing  method 
is  restricted  to  linear  systems.  Therefore,  we  apply  the  balancing  method  to  the  linearized  version  (56)  of 
the  nonlinear  wafer  heat  transfer  ODE  (19). 

Numerical  Difficulties  for  Nonminimal  Systems 


Numerical  difficulties  may  arise  in  the  use  of  the  balancing  model  reduction  procedure.  The  problem 
of  calculating  the  balancing  transformation  T^i  will  be  badly  conditioned  when  the  matrix  WCW0  has  a 
high  condition  number,  i.e.,  when  some  modes  are  nearly  uncontrollable  or  unobservable.  To  see  this, 
consider  Glover’s  algorithm  for  computing  the  balancing  transformation.  Given  Wc  and  WQ,  let  the  Cholesky 
factorization  of  WQ  be  given  by 

Wo  =  rtr. 

Difficulties  immediately  appear  when  W0  is  singular,  or  nearly  singular.  If  R  can  be  computed,  then  RWCRT 
is  diagonalized  as 

RWCR  =  UY?Ut 

with  UTU  =  1  and  S  =  diag[<Ji, . . .  <rn].  The  balancing  transformation  is  given  by 


Thai 


S  ~1/2UtR. 


Once  again,  difficulties  appear  when  Wc  is  singular,  or  nearly  singular. 

It  turns  out  that  we  encounter  precisely  these  difficulties  when  attempting  to  calculate  a  balancing 
transformation  for  the  linearized  wafer  temperature  field  evolution  system.  The  condition  numbers  of  the 
objects  of  interest  are 

cond(Wc)  =  3.2  x  1018 

condfff'o)  =  3.8  x  1018 

cond (WCW0)  =  9.4  x  1018 

so  that  the  system  is  nearly  uncontrollable  and  unobservable,  at  least  to  an  order  of  numerical  precision 
that  makes  it  impractical  to  compute  a  balancing  transformation  using  the  standard  technique.  This  result 
is  expected,  since  lamp  influence  functions  and  initial  wafer  temperature  profiles  are  always  smooth  and 
relatively  fiat.  Hence,  there  will  be  non-smooth  or  spatially  fluctuating  temperature  profiles  that  are  almost 
impossible  to  produce  using  the  available  control  inputs. 

Schur  Method  For  Balancing 
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In  order  to  alleviate  the  numerical  difficulties,  we  use  a  procedure  presented  by  Safanov  and  Chiang  [40] 
for  calculating  a  realization  of  a  /cth-order  reduced  model.  The  reduced  model  is  not  necessarily  balanced, 
but  has  transfer  function  G(s)  which  is  exactly  the  same  as  that  obtainable  by  applying  Moore’s  (or  Glover’s) 
balanced  realization  model  reduction  procedure  to  any  minimal  realization  of  G(s),  the  transfer  function  of 
the  original  nth-order  model.  Since  the  reduced  transfer  function  is  the  same,  it  enjoys  the  same  attractive 
error  bound  that  was  shown  by  Glover. 

The  procedure  presented  in  [40]  is  referred  to  by  Safanov  and  Chiang  as  a  “Schur  method”  for  balanced- 
truncation  model  reduction  because  the  resulting  transformation  is  computed  using  the  Schur  decomposition 
of  the  matrix  WCWQ .  The  idea  is  that  balancing  can  be  avoided  altogether.  Instead,  we  only  need  to 
compute  bases  for  the  left  and  right  eigenspaces  associated  with  the  “big”  eigenvalues  of  WCWQ.  The  ordered 
Schur  form  of  WCWQ  is  used  to  compute  orthonormal  bases  enabling  numerically  robust  computation  of  a 
nonbalanced  realization  of  the  /cth-order  reduced  model  G(s).  The  algorithm  is  stable  and  effective  without 
regard  to  nearness  to  unobservability  or  uncontrollability. 

General  Procedure 


The  general  procedure  is  summarized  as  follows.  Let  (A,  B ,  C,  D )  be  a  realization  of  G(s)  and  let  k 
denote  the  order  of  the  reduced  model. 


1.  Compute  matrices  Vr^ig^iMg  ^  Hnx/c  whose  columns  form  basis  for  the  respective  right  and  left 
eigenspaces  of  WCW0  associated  with  the  “big”  eigenvalues  g\ , . . . ,  o\.  This  is  done  using  the  ordered 
Schur  decomposition  of  WCWQ  and  is  described  in  the  next  paragraph. 


2.  Let 

and  compute  the  SVD 


Ebig  —  Vl,bigVr,big 
T 


Ue , big  £ E,big  ^E,big  • 

3.  The  not  necessarily  balancing  transformations  are 


SiMg  =  ViMgUBMa^tS,  e  RnXfc 

SrMa  =  VrMgVEMa^big  €  RnXfe 

so  that  the  reduced  model  is  given  by 

(A,  B,  C ,  D)  =  (S^bigASrMg,S^bigB,  CSrMg ,  D) 


Specific  Procedure  for  Eigenspaces  of  WCW0 


1.  Compute  an  orthogonal  real  matrix  V  such  that  VWCW0VT  is  upper  triangular,  i.e.,  put  WCWQ  in 
Schur  form.  The  fact  that  Wc  and  WQ  are  real  and  symmetric  ensures  the  existence  of  a  real  Schur 
transformation  matrix  V. 


2.  Compute  orthogonal  real  transformations  Va  and  Vd  which  order  the  Schur  forms  in  ascending  and 
descending  order,  respectively, 

V?WcW0Va  =  diag(Aa, A  a,i)  +  Ta  (59) 

VlWcWQVd  =  diag(Ad, Xd,n)  +  Td  (60) 

(61) 


where  Ta  and  Td  are  strictly  upper  triangular  and  such  that 


—  •  •  •  5  ^d,k}  —  {04  ,  •  •  •  , 

=  lj  •  •  •  5  ^d,n}  =  •  •  •  5  & n  } 


•  •  •  5  ^a,k  } 
{^a,/c+l5  •  •  •  5  ^a,n} 
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(62) 

(63) 

(64) 


Figure  40:  Left  basis  elements  for  balancing  transformation. 


Figure  41:  Right  basis  elements  for  balancing  transformation. 


3.  Partition  Va  and  Vd  as 


Va 


Vd 


(65) 

(66) 


Balancing  Results 

The  Schur  method  for  balancing  is  applied  to  the  linearized  model  (56)  and  (54),  i.e.,  the  realization 
(A,  B ,  C,  0).  The  condition  numbers  for  all  of  the  matrices  used  in  the  procedure  are  less  than  1000.  The  left 
and  right  basis  elements  for  the  balancing  transformation  are  shown  in  Figures  40  and  41.  The  corresponding 
relative  energy  values  are  given  in  Table  4. 

3.5  Comparison  and  Remarks 

Validation  of  the  predicitive  capability  of  the  reduced  models  is  accomplished  by  comparing  simulation 
results  using  the  original  nth-order  model  with  simulation  results  using  reduced  /cth-order  approximations 
for  various  values  of  the  reduced  model  order  k.  In  particular,  the  maximum  deviation  of  the  output  signals, 
i.e.,  the  thermocouple  readings,  between  the  original  and  reduced  models  are  computed  for  each  of  the  model 
reduction  approaches  we  have  previously  described. 
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Procedure 


Model  reduction  is  achieved  by  projecting  the  evolution  equations  of  our  wafer  thermal  model  onto  a  k- 
dimensional  subspace  (, k  <C  n)  of  the  original  n-dimensional  state  space.  This  is  done  via  standard  Galerkin 
projection.  The  ^-dimensional  subspace  is  chosen  as  the  span  of  a  set  of  basis  vectors  {0i, . . . ,  (/>&}  generated 
using  one  of  the  previously  described  POD  or  balancing  methods.  In  particular,  they  are  the  k  basis  elements 
with  the  largest  corresponding  energy  values.  Thus,  the  n  —  k  state  components  with  the  smallest  energy 
values  are  deleted  by  this  transformation.  Each  projection  results  in  a  /cth-order  ODE  model  whose  evolution 
in  time  yields  an  approximation  to  the  evolution  of  the  original  nth-order  model.  Finally,  an  approximation 
to  the  wafer  temperature  field  Tw  is  reconstructed  by  transforming  back  to  the  original  state  space. 

If  we  arrange  the  basis  vectors  in  a  matrix 

$  =  [4>  i . . .  oa'| 

then  the  approximation  to  the  wafer  temperature  field,  Tw,  is  given  by 

Tw  =  $z  (67) 

where  z  G  span{0i, . . . ,  (/>&}  is  the  state  variable  in  the  /cth-order  ODE  model 

z  =  $t[Ac$z  +  Ar(<f>z )4  +  Av<f>z  +  r  +  BP]  (68) 

resulting  from  Galerkin  projection.  Observe  that  <I>T  maps  elements  of  the  original  n-dimensional  state  space 
to  the  reduced  order  ^-dimensional  state  space,  while  <f>  maps  elements  of  the  reduced  order  space  back  to 
the  original  space.  Also  note  that  only  in  the  case  when  n  =  k  does  represent  a  complete  orthonormal 
basis,  i.e.,  <h<I>T  =  4>T4>  =  1.  Once  we  have  selected  4>,  using,  e.g.,  POD  or  balancing,  (68)  is  numerically 
integrated  and  an  approximation  to  Tw  is  reconstructed  using  (67). 

Validation  Tests 

The  model  reduction  methods  are  tested  by  numerically  integrating  the  original  model  and  reduced 
models  using  a  uniform  initial  temperature  field  and  two  test  recipes  as  lamp  control  inputs.  The  test 
recipes  are  different  from  those  used  to  generate  the  RSC  and  POC  data  ensembles. 

Test  Recipe  1  P  =  [0.5  0.5  0.5]  t  G  [0, 1] 

Test  Recipe  2  P  =  [1.0  0.0  0.0]  t  G  [0,0.4) 

P  =  [0.0  1.0  0.0]  t  G  [0.4,  0.7) 

P  =  [0.0  0.0  1.0]  t  G  [0.7, 1.0] 

The  reduced  models  for  k  =  1,  2,  3,  4,  and  5  are  numerically  integrated  using  the  same  control  recipes 
and  ICs  as  for  the  original  n  =  101  order  model.  Simulated  thermocouple  readings  are  recorded  for  each 
simulation.  The  error  between  the  original  and  reduced  models  is  computed  as 

e(fc)  =  ][ Ttc  -  Ttc||ma*  ,  k  =  1,  2, 3, 4,  5  (69) 

where  we  define  the  norm  ||i/||maa.  for  time-dependent  p- vector  y(t)  as 

IMImaz  =  max  {yi(t)  :  0  <  t  <  oo  ,  1  <  i  <  p]  (70) 

where  yi(t)  corresponds  to  the  temperature  reading  of  thermocouple  i  at  time  t.  Thus,  (69)  gives  the 
maximum  deviation  between  actual  and  estimated  thermocouple  readings  over  the  entire  simulated  time 
sequence  and  over  all  three  thermocouples,  i.e.,  a  “worst  case”  error. 

Results 

Due  to  the  shape  of  the  lamp  heat  flux  intensity  profiles  and  the  smoothing  effect  of  the  diffusion  operator, 
the  evolution  of  the  wafer  temperature  field  does  not  produce  especially  interesting  behavior,  e.g.,  spatial 
profiles  whose  fluctuations  from  the  mean  vary  substantially  in  the  mean  square  sense  from  the  initial  profile, 
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Percent  Energy  Associated  With  Transformation  Basis  Elements 


Method 

Mode  1 

Mode  2 

Mode  3 

Mode  4 

Mode  5 

POD  RSC 

95.06 

4.77 

0.14 

0.03 

0.00 

POD  POC 

93.43 

6.25 

0.23 

0.09 

0.01 

Balancing 

98.02 

1.83 

0.13 

0.02 

0.00 

Table  4:  Normalized  eigenvalues,  i.e.,  percent  energy,  corresponding  to  basis  elements  used  in  model  reduction 
for  POD  method  with  RSC  data,  POD  method  with  POC  data,  and  balancing  approach. 


Maximum  deviation  (degrees  C)  between  outputs  of  original  and  reduced  models 


Reduction 

Reduced  Model  Order 

Simulation 

Method 

1 

2 

3 

4 

5 

Test  1 

POD  RSC 

27.23 

2.68 

0.58 

0.11 

0.01 

POD  POC 

26.85 

1.26 

1.13 

0.10 

0.05 

Balancing 

50.68 

7.03 

0.44 

0.08 

0.02 

Test  2 

POD  RSC 

72.33 

5.22 

1.48 

0.18 

0.05 

POD  POC 

72.60 

4.79 

4.35 

0.43 

0.10 

Balancing 

80.81 

14.28 

1.70 

0.12 

0.04 

Table  5:  Maximum  deviation  (degrees  C)  between  outputs  of  original  and  reduced  models  for  POD  method 
with  RSC  data,  POD  method  with  POC  data,  and  balancing  approach. 


assuming  the  initial  profile  is  relatively  smooth.  Thus,  we  expect  little  difficulty  in  capturing  the  essence  of 
the  input-output  behavior  of  the  system  in  a  low  dimensional  model.  Our  results  show  that  this  is  indeed 
the  case. 

Tables  4  and  5  give  the  relative  energy  values  for  basis  elements,  and  the  maximum  thermocouple  tem¬ 
perature  deviations  for  the  original  and  reduced  order  models.  Figures  42,  43,  and  44  show  simulated 
thermocouple  readings  resulting  from  test  recipe  1,  for  the  original  n  =  101  order  model,  and  reduced  mod¬ 
els  of  order  k  =  1,  2,  and  3.  Figures  45,  46,  and  47  show  simulated  thermocouple  readings  resulting  from 
test  recipe  2.  From  the  figures,  it  is  clear  that  the  measured  response  for  the  reduced  models  is  close  to  that 
of  the  original.  We  further  analyze  the  data  as  follows. 

For  the  inputs  used  in  the  validation  tests,  the  input-output  behavior  of  the  wafer  heat  transfer  system 
can  be  reconstructed  using  reduced  models  of  order  4  so  that  thermocouple  readings  are  within  1  degree  C 
of  the  readings  using  the  original  model.  This  holds  whether  the  POD  or  balancing  method  is  used,  and  for 
whichever  set  of  empirical  data  was  used  for  computing  the  POD  transformation.  Even  reduced  models  of 
order  2  produce  a  reasonable  approximation  with  “worst  case”  errors  less  than  15  degrees  C. 

The  POD  method  appears  to  have  performed  slightly  better  than  the  balancing  method  in  this  study. 
One  reason  for  this  result  is  that  the  balancing  transformation  was  computed  for  the  linearized  system, 
while  the  validation  tests  were  performed  for  the  reduced  order  nonlinear  system.  Another  reason  is  the 
simple  input-output  behavior  of  this  particular  system.  The  principal  components  of  the  flow  are  relatively 
insensitive  to  the  choice  of  inputs,  and  hence,  any  set  of  principal  components  derived  from  empirical  data 
for  this  system  will  likely  be  efficient  for  model  reduction  purposes. 

3.6  Reduction  from  CFD  Models 

Figure  48  shows  time  series  data  points,  or  snapshots,  of  the  temperature  field  in  the  process  chamber  of 
the  ASM  Epsilon-1.  This  is  a  sampled  data  representation  of  the  gas  temperature  evolving  in  time  over  a 


49 


Original:  n  =  101 


Reduced: k  = 1 


Reduced:  k  =  2  Reduced:  k  =  3 


Figure  42:  Thermocouple  readings  for  original  and  reduced  models  with  Test  Recipe  1  using  transformation 
from  POD  RSC. 


Original:  n  =  1 01  Reduced:  k  =  1 


Figure  43:  Thermocouple  readings  for  original  and  reduced  models  with  Test  Recipe  1  using  transformation 
from  POD  POC. 


Original:  n  =  1 01  Reduced:  k  =  1 


Figure  44:  Thermocouple  readings  for  original  and  reduced  models  with  Test  Recipe  1  using  balancing 
transformation. 
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Original:  n  =  1 01  Reduced:  k  =  1 


Figure  45:  Thermocouple  readings  for  original  and  reduced  models  with  Test  Recipe  2  using  transformation 
from  POD  RSC. 


Original:  n  =  1 01  Reduced:  k  =  1 


Figure  46:  Thermocouple  readings  for  original  and  reduced  models  with  Test  Recipe  2  using  transformation 
from  POD  POC. 


Original:  n  =  1 01  Reduced:  k  =1 


Figure  47:  Thermocouple  readings  for  original  and  reduced  models  with  Test  Recipe  2  using  balancing 
transformation. 
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Time  =  0.02  seconds 


1000 
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Time  =  0.04  seconds 
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Time  =  0.06  seconds 
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1000 
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1000 
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Figure  48:  Snapshots  of  temperature  field  in  process  chamber:  t  =  0.02,  0.04,  0.06,  0.08,  0.10  seconds. 
Flowing  gases  are  heated  as  they  pass  the  wafer. 


period  of  0.1  seconds  under  typical  operating  conditions. 

In  order  to  generate  this  simulation  data,  the  fully  coupled  set  of  nonlinear  transport  equations  were 
solved  on  a  grid  with  720  nodes  using  Fluent.  Even  if  the  energy  balance  equations  for  gas  temperature 
were  uncoupled  from  the  fluid  flow  and  species  transport  equations,  the  simulator  is  numerically  integrating 
a  model  with  720  states,  one  for  each  grid  point.  Coupling  among  transport  equations  makes  the  order  of 
the  system  much  higher  than  that. 

This  simulation  required  approximately  15  seconds  of  computing  time  to  converge  to  a  solution  for  each 
simulated  time  step,  or  a  factor  of  750  slower  than  the  actual  time  period  being  simulated.  This  means 
that  not  only  is  the  CFD  model  computationally  demanding  for  simulation,  but  also  it  is  computationally 
prohibitive  for  purposes  of  real-time  model-based  control.  In  addition,  controller  design  and  optimization 
for  a  system  with  hundreds  of  states  will  be  excessively  complicated. 

We  apply  our  methodology  for  model  reduction,  in  which  we  extract  spatial  structures,  or  patterns,  that 
dominate  the  flow  dynamics  as  given  by  the  snapshots  of  the  heat  flow.  These  “coherent  structures” ,  or  prin¬ 
cipal  components,  are  shown  in  Figure  49.  We  note  that  these  are  the  empirically  determined  eigenfunctions 
of  the  two-point  spatial  covariance  of  the  fluctuations  from  the  mean  heat  flow.  They  contain  92.7%,  6.9%, 
and  0.5%,  respectively,  of  the  energy  in  this  flow.  This  preliminary  analysis  suggests  that  we  can  model  the 
evolution  of  gas  temperature  with  a  model  of  dimension  2  to  a  high  degree  of  correctness. 

A  similar  analysis  has  been  done  with  time  evolution  of  silane  mole  fraction.  Figure  50  shows  the  snapshots 
recorded  over  a  simulated  period  of  0.2  seconds  (silane  mole  fraction  reaches  steady-state  in  approximately 
double  the  time  of  gas  temperature).  Figure  51  shows  the  corresponding  principal  components.  Again,  99% 
of  the  energy  is  contained  in  2  basis  functions. 

4  Conclusion 

Development  of  high-fidelity  and  low-order  models  is  the  first  step  in  achieving  the  objectives  of  this  project. 
The  task  now  is  to  use  the  results  to  improve  manufacturing  effectiveness.  To  this  end,  we  outline  here  some 
directions  for  further  work,  and  summarize  what  we  have  accomplished  so  far. 
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Energy  =  92.67% 


Energy  =  6.86% 


Figure  49:  Principal  components  (empirically  determined  eigenfunctions)  of  heat  flow  in  process  chamber, 
corresponding  to  energies  92.7%,  6.9%,  and  0.5%  (top  to  bottom). 


Time  =  0.04  seconds 


Time  =  0.08  seconds 


Time  =  0.20  seconds 


0.02 

0.015 

0.01 

0.005 


0.02 

0.015 

0.01 

0.005 


0.02 

0.015 

0.01 

0.005 


0.02 

0.015 

0.01 

0.005 


0.02 

0.015 

0.01 

0.005 


Figure  50:  Snapshots  of  silane  mole  fraction  in  process  chamber:  t  =  0.04,  0.08,  0.12,  0.16,  0.20  seconds. 
Silane  is  depleted  as  surface  reactions  take  place. 
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Energy  =  94.68% 


Energy  =  4.77% 


Figure  51:  Principal  components  (empirically  determined  eigenfunctions)  of  silane  species  transport  in  pro¬ 
cess  chamber,  corresponding  to  energies  94.7%,  4.8%,  and  0.5%  (top  to  bottom). 

4.1  Future  Directions 

There  are  several  ways  in  which  we  can  attempt  to  improve  manufacturing  effectiveness.  One  is  by  improving 
product  quality  as  measured  by  deposition  uniformity.  Another  is  by  providing  the  process  engineer  with 
tools  to  increase  efficiency  and  provide  insight  and  predictions. 

Of  particular  interest  is  the  improvement  of  deposition  uniformity  in  the  presence  of  microfeatures  on 
the  wafer  surface.  One  way  in  which  we  hope  to  attack  this  problem  is  by  quantifying  the  tradeoffs  between 
pattern  pitch  and  uniformity.  Then,  for  a  given  uniformity  requirement,  we  can  find  the  minimum  allowable 
pattern  pitch,  or  for  a  given  pattern  pitch,  we  can  find  the  maximum  possible  uniformity. 

Currently,  the  mechanisms  which  affect  deposition  uniformity  in  the  presence  of  oxide  patterns  on  the 
wafer  surface  are  not  well  understood.  Development  of  the  feature  scale  model,  which  will  work  in  conjunc¬ 
tion  with  the  reactor  and  process  model,  will  provide  insight  toward  understanding  these  mechanisms  and 
important  cause  and  effect  relationships. 

A  process  recipe  consists  of  a  sequence  of  processing  steps,  each  of  which  performs  a  given  task  over  a 
fixed  time  period,  and  is  characterized  by  a  set  of  operating  conditions  selected  by  the  process  engineer.  Such 
operating  conditions  include  choice  of  process  gases,  inlet  flow  rates,  chamber  pressure,  and  wafer  temper¬ 
ature.  The  ASM  Epsilon-1  reactor  system  has  the  capability  of  controlling  these  conditions  automatically, 
according  to  the  pre-programmed  recipe,  during  the  processing  run.  In  addition,  there  are  several  settings 
over  which  the  process  engineer  has  control  and  which  remain  constant  throughout  the  process  recipe.  These 
are  manually  set  by  the  engineer,  and  include  inlet  injector  slit  widths,  thermocouple  offsets,  and  PID  control 
gains  in  temperature  control  loops. 

A  given  process  recipe  and  set  of  reactor  settings  will  produce  a  corresponding  thin  film  deposition  on 
the  wafer  surface.  We  are  interested  in  certain  characteristics  of  this  result  including  deposition  thickness, 
uniformity,  and  morphology.  Therefore,  we  want  our  models  to  ultimately  have  the  capability  to  predict 
these  characteristics  given  a  process  recipe  and  set  of  reactor  settings.  Once  this  capability  is  achieved, 
and  experimentally  verified,  the  models  can  be  used  for  development  of  tools  to  aid  the  process  engineer  in 
determining 

•  optimal  thermocouple  offsets  for  a  given  uniformity  requirement  and  a  given  set  of  process  conditions; 
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•  optimal  inlet  injector  slit  widths  for  a  given  uniformity  requirement  and  a  given  set  of  process  conditions; 

•  final  thickness,  deposition  rate,  film  composition,  uniformity,  and  morphology  for  a  given  process  recipe 
and  given  thermocouple  offsets  and  inlet  injector  widths,  all  relevant  recipe  steps  included;  input  to 
this  tool  would  be  similar  to  how  the  recipe  is  input  to  the  ASM  Epsilon- 1  itself;  output  would  be  a 
set  of  predicted  results; 

•  optimal  pressure,  flow  rate,  and  temperature  recipes; 

•  deliterious  effects  such  as  deposition  on  walls  and  mechanical  stresses  on  equipment  and  product  for  a 
given  recipe; 

•  optimal  PID  gains  for  Foxboro  temperature  controllers. 

The  tools  would  be  useful  in  ”what  if’  analysis  of  various  production  runs.  The  recipes  could  be  tested 
using  simulation  without  actually  performing  the  run,  avoiding  waste  of  time,  energy,  and  materials,  and 
allowing  for  rapid  tuning  of  adjustable  parameters  and  testing  of  ideas  before  production.  The  process 
engineer  could  verify  that  the  thermocouple  offsets  and  inlet  slit  widths  being  used  are  good  or  find  ones 
that  provide  better  uniformity.  The  tools  should  also  provide  some  insight  as  to  why  certain  settings  are 
better  than  others. 

Successful  implementation  of  low-order  process  and  equipment  models  will  also  be  useful  in  design  of 
real-time  model-based  control  strategies  for  tracking  and  regulation  of  desired  temperature  profiles.  Using 
sophisticated  control  systems  simulation  software  packages  such  as  Simulink4  and  SystemBuild5,  “system 
blocks”  can  be  constructed  from  the  low-order  models.  These  system  blocks  interface  to  other  system  blocks 
via  user-defined  inputs  and  outputs.  In  this  manner  a  controller  block  can  be  designed  and  a  feedback 
control  strategy  can  be  simulated.  The  power  of  these  models  can  be  extended  further  via  commercially 
available  rapid  prototyping  tools,  which  work  in  conjunction  with  the  simulation  software,  and  provide  the 
capability  for  system  blocks  to  connect  physically  with  the  actual  inputs  and  outputs  of  the  systems  they 
are  simulating. 

4.2  Summary 

We  have  described  work  done  to  date  on  a  joint  project  between  Northrop  Grumman  ESSD  and  the  ISR 
to  investigate  the  epitaxial  growth  of  Si-Ge  heterostructures  in  the  ASM  Epsilon-1  RTCVD  reactor.  An 
important  goal  of  this  research  is  to  improve  manufacturing  effectiveness  by  controlling  and  optimizing 
deposition  characteristics  such  as  uniformity  in  the  presence  of  microfeatnres  on  the  wafer  surface. 

High-fidelity  models  have  been  developed  and  implemented  using  sophisticated  general  purpose  and 
specialized  CFD  software  as  well  as  stand-alone  hard-coded  programs.  The  models  describe  transport 
mechanisms  and  chemical  reactions  in  the  reactor  including  heat  transfer,  fluid  flow,  species  transport,  and 
deposition  kinetics.  Experimental  work  has  been  done  to  determine  growth  parameters  and  validate  lamp 
heating  models. 

A  comparison  of  methods  for  model  reduction  via  linear  transformation  and  truncation  has  been  per¬ 
formed.  Low  order  models  derived  using  both  the  POD  and  balancing  methods  appear  to  provide  an  accurate 
approximation  to  the  original  model.  The  input-output  behavior  of  the  heat  transfer  system  can  be  recon¬ 
structed  to  a  high  degree  of  correctness  with  models  of  dimension  4  and  higher.  The  wafer  heat  transfer 
model’s  lack  of  “interesting”  behavior  facilitates  the  model  reduction,  and  causes  there  to  be  little  difference 
between  approximation  properties  of  reduced  models  derived  from  the  POD  and  balancing  approaches  in 
this  case.  Preliminary  results  in  applying  model  reduction  techniques  to  other  mechanisms  such  as  transport 
of  energy  in  the  gas  phase  and  transport  of  silane  species  suggest  that  low-order  models  can  be  developed 
to  accurately  predict  their  time  evolution. 

Efforts  toward  incorporation  of  more  complicated  mechanisms,  integration  of  the  various  models,  and 
model  reduction  will  lead  to  the  development  of  useful  tools  for  improvement  of  manufacturing  effectiveness 
and  product  quality. 


4The  Mathworks,  Inc.,  Natick,  MA 

5 Integrated  Systems,  Inc.,  Sunnyvale,  CA 
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Appendices 


A  Physical  Constants 

Listed  here  are  the  physical  constants  used  in  the  models.  The  units  have  been  selected  for  convenience  and 
consistency.  Properties  of  the  wafer  are  those  of  pure  silicon.  Chamber  wall  properties  are  those  of  quartz. 
Properties  of  the  process  gases  are  those  of  hydrogen  at  1000  K  and  1  ATM.  Chemical  kinetics  parameters 
are  those  experimentally  determined  from  reactions  involving  thermally  activated  deposition  of  polysilicon 
from  30  seem  of  2%  silane  in  hydrogen. 


Constant 

Description 

Value 

Units 

k0 

Arrhenius  Coefficient 

3.0787  x  103 

cm  sec-1 

Ea 

Activation  Energy 

1.6330  x  105 

J  mol-1 

Rg 

Gas  Constant 

8.314 

J  mol"1  K"1 

href 

Reference  Thickness 

1.0  x  10“4 

cm 

fir 

Rate  Pre-Exponential  Constant 

1.8472  x  10y 

dimensionless 

fie 

Rate  Exponential  Constant 

2.8059  x  101 

dimensionless 

kw 

Thermal  Conductivity  of  Wafer 

0.22 

W  cm-1  K-1 

Pw 

Mass  Density  of  Wafer 

2.3 

g  cm-3 

Cpw 

Heat  Capacity  of  Wafer 

2.3 

J  g"1  K-1 

Vb 

Boltzmann  Constant 

5.677  x  W71 

W  cm-2  K-4 

Emissivity  of  Wafer 

0.7 

dimensionless 

Qiw 

Absorptivity  of  Wafer 

0.5 

dimensionless 

Rw 

Radius  of  Wafer 

7.62 

cm 

A, 

Thickness  of  Wafer 

0.05 

cm 

hv 

Convective  Heat  Transfer  Coefficient 

2.6474  x  10“4 

W  cm-2  K-1 

Re 

Reynolds  Number  of  Gas  Flow 

27.2 

dimensionless 

kg 

Thermal  Conductivity  of  Gas 

4.40  x  10“3 

W  cm-1  K-1 

Pr 

Prandtl  Number  of  Gas  Flow 

0.686 

dimensionless 

L 

Chamber  Length 

50.8 

cm 

T 

-L  C 

Chamber  Wall  (Ambient)  Temperature 

700 

K 

ec 

Emissivity  of  Chamber  Wall 

0.37 

dimensionless 

T 

±9 

Gas  Temperature 

300 

K 

T 

Reference  Time 

60 

seconds 

Qref 

Reference  Heat  Flux 

29.24 

W  cm-2 

B  Dependent  Variables 


Variable 

Parameters 

Description 

Units 

Dimensionless  Conversion 

T 

-1-  W 

t,r 

Wafer  Temperature  Field 

K 

TW  /  Tamb 

XsiH4 

r 

Silane  Mole  Fraction 

dimensionless 

XsiHA 

CsiH4 

r 

Silane  Concentration 

mol  cm-3 

CsiH4  /  Ctot 

PsiHA 

r 

Silane  Partial  Pressure 

torr 

PsiHA  /  Ptot 

Rsi 

t,r 

Silicon  Deposition  Rate 

cm  sec-1 

h 

t,r 

Film  Thickness 

cm 

h  j 

Qi 

r 

Lamp  i  Heat  Flux  Profile 

W  cm-2 

Qi  /  Qref 

Ui 

t 

Lamp  i  Power  Setting 

dimensionless 
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C  Independent  Parameters 


Parameter 

Description 

Units 

Domain 

Dimensionless  Conversion 

t 

Time 

seconds 

[0, 00 ) 

f /r 

r 

Radial  Position 

cm 

[0,-Rto] 

r  /  Rw 

e 

Azimuthal  Position  (Angle) 

radians 

[0,2tt] 

9  /  2tt 

z 

Axial  Position  (Height) 

cm 

[0,A*] 

H 

< 

D  Balancing  For  Linear  Systems 

Consider  the  nth-order  stable  linear  system  with  minimal  state  space  realization  (A,  B,C,D)  of  transfer 
function  G{s)  =  D  +  C(s'SL  —  A)~1B  and  associated  controllability  and  observability  Gramian  matrices,  Wc 
and  W0  respectively,  determined  by  the  Lyapunov  equations 

WCAT  +  AWC  +  BBt  =  0  (71) 

and 

W0A  +  AtWc  +  CTC  =  0.  (72) 

It  is  known  [33,  12]  that  there  exists  an  invertible  linear  transformation  of  the  state  space  Thai  G  R"  x "  such 
that  the  transformed  realization 

(Abal ,  Bfral ,  Cbal ,  Dbal )  =  ATbal ,  B ,  CTbal ,  D) 

has  controllability  and  observability  Gramian  matrices  of  the  form 

Wc-bai  =  ^1Wc(rb-i1)T  =  diag(S1,S2,0,0) 

Wo-bai  =  T^alW0Tbai  =  diag(Si,  0,  S3, 0) 

where  Ei,  £2,  ^3  are  positive  definite  diagonal  matrices  given  by 

Si  =  diag(<r  1 ,  ...,ar) 

O 1  ^  O 2  O  >p  ^  <7r_|_l  —  ’  ’  ’  —  O yi  —  0 

and  the  cq  are  determined  by 

^  =  (A  i(WcW0))1/2. 

Such  a  transformation  is  called  a  balancing  transformation  and  the  resulting  state  space  realization  is  called  a 
balanced  realization.  When  transformed  to  the  balanced  realization,  the  input-output  influence  of  each  state 
component  relative  to  the  others  is  indicated  by  the  corresponding  singular  value  cq.  Thus,  the  realization  can 
be  partitioned,  and  the  n  —  k  least  reachable-observable  states  truncated,  yielding  the  /Th-order  transfer 
function  G(s).  It  has  been  shown  that  the  reduced  model  in  stable,  minimal,  and  balanced,  i.e.,  both 
controllability  and  observability  Gramian  matrices  are  equal  to 

£ bal  =  diag((7 1 ,  .  .  .  ,crfc). 

It  has  also  been  shown  that  the  reduced  model  enjoys  the  error  bound 

r 

o(G(j(jo)  —  G{j(jj))  <  2  ^  Oi  \/u.  (75) 

i=k+ 1 


(73) 

(74) 
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